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SUMMARY 


An inverse swept wing code is described that is based on 
the widely used transonic flow program FL022. The new code 
incorporates a free boundary algorithm permitting the pressure 
distribution to be prescribed over a portion of the wing surface. 
A special routine is included to calculate the wave drag, which 
can be minimized in its dependence on the pressure distribution. 
An alternate formulation of the boundary condition at infinity 
has been introduced to enhance the speed and accuracy of the 
code. A FORTRAN listing of the code and a listing of a sample 
run are presented. There is also a user's manual as well as 
glossaries of input and output parameters. 


INTRODUCTION 


After much controversy about shockless airfoils in the 
theory of transonic flow, experimental work has established that 
wings can be constructed to virtually eliminate wave drag over a 
practical range of supercritical speeds. Computational fluid 
dynamics has become a primary tool for the design and analysis 
of these supercritical wings. More specifically, computer codes 
developed at NYU to calculate transonic flow in both two and 
three dimensions have become widely accepted by the aircraft 
industry. It is our purpose here to describe and list the latest 
of these codes, which serves to redesign a swept wing by select- 
ing its pressure distribution so that the wave drag is minimized 
at a fixed speed and angle of attack. 

Perhaps the best way to design shockless airfoils in two- 
dimensional transonic flow is to use the hodograph transformation 
in combination with analytic continuation into the complex 
domain [1,3]. Most analysis codes, on the other hand, depend on 
an introduction of artificial viscosity and artificial time that 
is motivated by the retarded difference scheme of Murman and Cole 
[9] . The method of artificial viscosity can also be applied to 
the design problem, for which it is especially helpful in three 
dimensions [6,7]. An approach of this kind has been adopted in 
developing the inverse swept wing code we are now concerned with. 
Our procedure has been to modify the FL022 code of Jameson and 
Caughey, which is in turn based on an earlier oblique wing code 
for the calculation of transonic flow in three dimensions [2,8]. 

In the next section of the report we shall review theoreti- 
cal aspects of the transonic s codes that are either somewhat 
controversial or have not been well publicized elsewhere. 
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Indications will be given of how the basic method can be general- 
ized; but detailed treatments of more complicated problems, such 
as the three-dimensional flow through a cascade of compressor 
blades, will be left to other publications. An example of a 
supercritical swept wing that has been redesigned by applying 
the new three-dimensional code will be discussed. Both computa- 
tional and physical properties of the example will be emphasized. 
Then a detailed description of the code will be presented that 
can serve as a user's manual. The final sections of the report 
are devoted to the listing of a sample run for the supercritical 
wing just referred to and to a listing of the code, with comment 
cards . 


MATHEMATICAL BACKGROUND 


The transonic flow around airfoils and wings is usually 
calculated by considering a velocity potential <j> that satisfies 
the second order quasilinear partial differential equation 


(c^-u^)^ + (c^-v^)<t> + (c^-w^)4> - 2uv(|> - 2vw(J) - 2wutj) = 0 , 

XX yy zz xy yz zx 

where u = ^ # v = <j) and w = <j) are the velocity components 

and c is the^speed ^ of sound ^defined by Bernoulli's law 


2 ^ 2 ^ 2 
U + V + w 



const . 


A Neumann problem for <}) is specified by setting its normal deriv- 
ative equal to zero at. the wing and prescribing its asymptotic 
behavior at infinity. Finite difference schemes that capture 
weak shock waves effectively are arrived at by adding an artifi- 
cial viscosity term to the equation for <i> . This term is obtained 
by retarding difference approximations to second derivatives in 
the direction of the flow, which does not alter the boundary 
condition at the wing. Iterative methods to solve the difference 
equations for <|) are found by introducing artificially time- 
dependent terms that force decay to a steady state [s]. 

An objection can be made to use of the velocity potential 
because that presumes constant entropy, whereas the wave drag, 
which is of primary interest, has the same order of magnitude 
as the jump in entropy across shocks, to which it can even be 
attributed. However, we have been able to develop an expression 
for the wave drag in terms of the velocity potential that is 
accurate to lowest order for weak normal shock waves [6] . This 
is important because there are ambiguities in determining a 
steady state solution of the Euler equations that are perhaps 
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best overcome by the assumption of irrotationality that charac- 
terizes potential flow. More general steady solutions may 
include vortices such as those that occur in models of the wake. 
Therefore some hypothesis must be made to ensure uniqueness of 
the flow. 

In two dimensions another possibility for handling the 
Euler equations is to introduce the stream function ^ as an 
independent variable and to calculate the flow by solving a 
partial differential equation for the ordinate y as a function 


y = y(x,<j^) 


of X and . This is equivalent to making the topological assump- 
tion that each vertical line intersects each streamline just 
once, which does eliminate vortices. But it is awkward to 
formulate the laws of conservation of momentum in a fashion 
convenient for numerical computation of the unknown y. Further- 
more, experience with analogous problems in the calculations of 
magnetohydrodynamic equilibrium shows that existence as well 
as uniqueness becomes questionable for steady solutions of the 
Euler equations in three dimensions. The analogy is based on 
letting the velocity, the vorticity and the Bernoulli constant 
for incompressible flow correspond respectively to the magnetic 
field, the current density and the pressure in magnetohydro- 
dynamics [5] . 

How the wave drag may be represented in terms of the velocity 
potential (p is most readily understood by studying a model problem 
for one-dimensional flow. Application of the retarded difference 
scheme of Murman and Cole to the small disturbance equation for (fi 
leads us to consider the ordinary differential equation 


describing conservation of mass, where h is a positive mesh size 
parameter multiplying a term on the right that we conceive of as 
artificial viscosity. The flow is said to be supersonic when 
c|)x > 0 and subsonic when <j) < 0 . If appropriate boundary condi- 

tions of the form 


<j)(a) = = C > 0 

are imposed at the ends of the interval [a,b] , a unique solution 
is found that approaches a pair of intersecting lines with the 
opposite slopes C and -C as h ^ 0 . The intersection of the lines 
is a shock wave across which remains continuous [2] . 
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Multiplying by <|)jj on both sides of the ordinary differen 
tial equation for (f> , we obtain an analogue 
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h[max(4>x/0)<!)^<|)^^]^ 


h 


max 



0 )(|) 


2 

XX 


of the law of conservation of momentum. Integration by parts 
and passage to the limit as h 0 yields the entropy inequality 

- I 

b 

= - lim h I max((|) dx < 0 . 

h-^o J ^ 

a 

This not only establishes the necessity of the requirement C > 0 
in our model problem, but also suggests that the integral on the 
right is a legitimate measure of both the wave drag and the jump 
in entropy. A similar argument has been used to represent the 
wave drag as a volume integral involving <j> for potential flow in 
both two and three dimensions [6,7]. The resulting formula has 
been implemented in our swept wing code and enables us to plot 
shock waves in a fashion indicating the amount of drag associated 
with them. 

It is important to realize that the integrand in the volume 
integral for the wave drag depends in a subtle way on the form 
of the artificial viscosity used to calculate (p . To understand 
why this should be so one has only to alter the artificial 
viscosity on the right in the ordinary differential equation 
given above for <P to obtain 


(|) (j) = h (|) 

X XX XXX 


instead. The same solution as before is found in the limit as 
h -»■ 0. The resulting entropy inequality 

b 

- $ = - lim h I dx < 0 

3 h»0 J 

a 

remains unaltered except that there is a change in the integrand 
on the right. Thus the integral representing the wave drag is 
seen to have the same value it had before, but the way in which 
the shock wave is smeared when h > 0 becomes significantly difi^ 
ferent. 

Another issue that arises in the computation of transonic 
flow around airfoils and wings is whether or not to put the 
finite difference equations in conservation form. Strictly 
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speaking this must be done to approximate the shock conditions 
accurately. However, boundary layer-shock wave interaction 
and, more specifically, the pressure recovery at the foot of a 
normal shock wave are poorly modeled by the conservation form 
of the equation for <j) . This seems to be due to a term in the 
artificial viscosity that can be eliminated by reverting to a 
simpler difference scheme that is closely related to the 
original method of Murman and Cole. We have chosen to retain 
such a nonconservative scheme in the swept wing code listed in 
this report. However, it is not difficult to modify the code 
so as to bring the equation for <j) into the mathematically more 
correct conservation form. 

For the model problem of one-dimensional flow the artificial 
viscosity in conservation form is 


h[4><|) ] =h<|)(() +h(|) 

X XX X X XXX XX 


whereas the nonconservative version is just h ‘('x'f^xxx • 
difference between these two viscosities is a positive term h<J)xx 
referred to above that represents mass generated by shock waves. 
For the full transonic flow problem the analogous quantity may 
contribute significantly to truncation error in supersonic 
regions where no shocks occur. Its omission from the nonconserva- 
tive scheme adopted in the swept wing code therefore has 
the advantage of improving accuracy to a certain extent on the 
crude meshes that one must resort to for a three-dimensional 
calculation of this kind. Moreover, the nonconservative scheme 
seems especially appropriate for flows that are designed to be 
shockless anyway. 

Our principal concern is the inverse problem of shaping a 
swept wing so that its pressure distribution may be prescribed. 
More specifically, we wish to choose the surface y = f(x,z) of 
the wing so that the square of the speed q assumes given values 
qo(x,z)2. This requirement yields a free boundary condition 

Q(f,f^,f^) = qQ(x,z)^ - q^ =0 

which we may view as a partial differential equation of the first 
order for the unknown function f. In the implementation of the 
computer code x, y and z are taken to be sheared parabolic 
coordinates such that the surface of the wing lies near the 
plane y = 0 and the flow is restricted to the half-plane y > 0. 
Problems with closure are circumvented by introducing a fixed 
surface y = fg(x,z) and imposing the constraiht f ^ fg , which 
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asserts that the wing must enclose a specified inner structure. 
The free boundary condition is only supposed to be fulfilled 
at points where f > fQ . Difficulties in locating stagnation 
at the leading edge or with closure at the trailing edge are 
avoided by choosing the assigned speed qQ so that it decays 
rapidly there and makes f = fg outside a range in the middle 
of the wing where the free boundary condition becomes operative. 

The free boundary problem we have formulated seems to be 
well posed even in the case of transonic flov/, but hanging 
shocks tend to appear above the wing even when the prescribed 
pressure distribution is smooth at the surface. An iterative 
scheme to solve the free boundary problem numerically is arrived 
at by letting the free surface function f vary suitably with 
the artificial time parameter t of the transonic flow calcula- 
tion. The motion of the free surface is defined by requiring 
that a partial differential equation of the fo 2 nn 


®l^xt ^2^t ^ 

be satisfied at points where f > fQ . The coefficients a- are 
adjusted to achieve convergence of the method. An analogue of 
the Lax-Wendroff finite difference scheme is used for the 
computation of f. Precisely how this is accomplished is a 
more subtle matter concerning which the reader is referred to 
the listing of the code. 

To eliminate shock waves on a swept wing at given speed 
and angle of attack it does not suffice just to prescribe the 
pressure smoothly. Experience with shockless airfoils designed 
by the hodograph method suggests that to suppress shocks the 
pressure distribution ought to be peaky at the front of the 
supersonic zone. We have incorporated in the code an exponen- 
tial spline routine that generates desirable distributions of 
speed depending on relatively few free parameters [7] . There 
is also an option that enables one to choose the prescribed 
distribution so as to minimize the wave drag. The code can be 
used to design swept wings that achieve virtually shockless 
transonic flow at a specified condition. While this requires 
some skill because the problem of design is far from easy in 
practice, the code does seem to be robust and is capable of 
delivering meaningful results for a relatively modest expenditure 
of computational resources . 

The speed and performance of the code have been improved 
by vectorization and by modifying the boundary condition at 
infinity. The new boundary condition is imposed on a control 
surface at some distance from the wing. It asserts that the 
reduced potential should decay at a specified rate as its 
argument approaches infinity. This is equivalent to a linear 
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combination of homogeneous Dirichlet and Neumann conditions on 
the control surface. The precise rate of decay has been 
adjusted semi-empirically to give optimal results. 

The computational methods used to develop the inverse 
swept wing code can be applied to a variety of harder problems. 
Of particular interest are the flow past propellers or through 
cascades of compressor bladds and the flow around airplane 
configurations that include a fuselage or engine nacelles. 

The most crucial issue is how to treat complicated geometries 
in space of three dimensions. While the question of adequate 
coordinate generation remains a challenge, it would appear that 
the difficulties associated with transonic flow are now well 
in hand. It may also be worth inquiring whether less directly 
related problems, such as that of ship wave resistance, might 
be attacked successfully by similar techniques of computational 
fluid dynamics. 


DESIGN OF A STANDARD SUPERCRITICAL WING 


Codes for the analysis of transonic flow around given 
bodies have been used quite successfully to simulate wind 
tunnel measurements, especially in the case of two-dimensional 
motion. The agreement between computed and experimental data 
is usually excellent, and the calculated results are obtained 
more quickly and sometimes more cheaply. In Fig. 1 we display 
a typical comparison of theoretical estimates of the drag 
coefficient Cq with experimental measurements that shows how 
well the new formula for the wave drag we have described works 
in practice for two-dimensional flow. 

Design codes have been received less enthusiastically by 
the engineering community. They leave more choices up to the 
user, and the outcome of the computations may be less tangible. 

In this section we present an example of a supercritical wing 
that has been redesigned through an application of the inverse 
swept wing code. The results serve primarily as a sample run 
to illustrate how the code works, but they are also not without 
physical interest despite the crudeness of the mesh that is used. 

It is best to construct the swept wing from a supercritical 
airfoil to begin with. For this purpose we have chosen a 13% 
thick airfoil that is shockless at free stream Mach number 
M = 0.75 and lift coefficient = 0.54. The wing is swept 
back through an angle of 30° and has aspect ratio A = 3.8. 

To redesign the wing, which has noticeable shocks near the 
wall, a typically shockless pressure distribution has been 
prescribed. It is specified in vertical sections by exponential 
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splines defined over three adjacent intervals. Optimization 
was used to select peak values of the pressure so that the 
wave drag became a minimum at the design condition of free 
stream Mach niimber M = 0.83 and lift coefficient = 0.41. 

This reduced the wave drag coefficient from 0.0028 to 
0.0008 and softened the shock waves perceptibly. 

Fig. 2 shows how the plot routine of the code displays 
the results of an analysis calculation for our wing, and Fig. 3 
shows corresponding data after the design run. Five sections 
of the wing are seen on the right, and corresponding distribu- 
tions of the pressure coefficient Cp over the upper surface 
are seen on the left. Shock waves are plotted above the wing 
with a thickness indicative of the wave drag associated with 
them. The detailed input and output of the design run are 
listed in the report. The mesh consisted of 128 10 12 points, 

and 100 iterations were performed to achieve acceptable conver- 
gence. This took 10 minutes of machine time on the CDC 6600 
computer. 


DESCRIPTION OF THE CODE 


The NYU inverse swept wing code can be used for both 
analysis and design of a swept wing. In the design mode an 
option is available which invokes an optimizer (Harwell Mathe- 
matical Subroutine Library, AERE, England) to minimize the wave 
drag. 


The analysis mode is like FL022, which calculates the 
transonic flow past a given swept wing [8]. For analysis, data 
for the code is input on cards and stored on Tape 5 or read 
directly into Tape 5. The input consists of computation param- 
eters, wing geometry, and physical specification of the flow. 
The resulting information is stored on Tape 7. This file can be 
saved and is used to initialize the computation if a wing is 
to be designed. 

In the analysis mode the principal difference between the 
NYU inverse swept wing code and FL022 is the introduction of a 
new boundary condition at infinity for the reduced potential G. 
This condition is imposed on an outer control surface and it 
improves the speed and accuracy of the computation. It has 
the form 


Gj^l = (1 - 3) G. , 

where the index j+1 refers to a ghost point just outside the 
computational domain. The positive parameter 3 is scaled so 
that the requirement models a mixed Dirichlet and Neumann 


8 



condition 


3G 


+ 



= 0 


that serves to annihilate the leading term 1/r of a hypothetical 
expansion of G in spherical harmonics. In practice 3 has been 
chosen semi- empirically to give optimal resolution. The new 
boundary condition provides a more effective way of asserting 
that G decays at infinity [4,7] . 

In the design mode the surface of a given wing is modified 
so that a prescribed Mach number distribution is assumed over a 
portion of the wing. The optimization package attempts to 
minimize the wave drag by changing the Mach number distribution 
systematically. The wave drag is computed using a formula that 
has been discussed in the section on mathematical background. 

More precisely, the factor 4>x occurring there is replaced 
by a term M^ - 1 involving the Mach number M, and the derivative 
4>xx is replaced by a second derivative <j>ss in the direction of 
the flow. This results in an expression of the form 

C.^TT = I max(M^-l,0) <j>^ 

DW ^ ss 

for the wave drag coefficient , where the summation is extended 
over all mesh points and A is a factor determined by the finite 
difference equations that are used. Contributions from rarefac- 
tion waves are automatically excluded by the code. 

To modify a given wing in the design mode, an analysis run 
for the unmodified, or baseline, wing is made to assess the 
characteristics of the flow field. The resulting speed distri- 
bution is examined and used to construct a more desirable distri- 
bution. The new distribution is input to the code on Tape 9. 

An exponential spline routine in the code calculates the desired 
distribution based on the input. This should have approximately 
the same spanwise distribution of lift as the original wing. 

It should be smooth throughout the supersonic zone, but may be 
peaky near the leading edge. In addition to the design distri- 
bution, a wing surface is prescribed that is identical to the 
original baseline wing near the leading and trailing edges but 
may be slightly thinner near the middle of each spanwise chord. 

The design scheme adds material to this underlying shape to fill 
in the thinned areas in a manner consistent with the assigned 
speed distribution. The thinning is done by introducing a groove. 
The parameters defining the shape of the groove are input to 
the code on Tape 9. To avoid difficulties with trailing edge 
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closure and maintenance of thickness-to-chord ratio, the surface 
modifications are made on a region of the upper surface that 
excludes the leading and trailing edges. 

The computational domain has been obtained by applying a 
square root transformation to the physical domain that results 
in a representation Y = SQ(X,Z) of the wing as a shallow bump 
lying above the (X,Z) -plane. The wing surface is changed 
iteratively, starting from the original shape as an initial 
guess. At each step one or more cycles of line relaxation are 
done in the standard analysis mode. The resulting speed distri- 
bution is compared with the desired distribution. Surface 
modifications are made that depend on the difference between 
the two distributions. If the modifications cause the computed 
surface to fall below the prescribed underlying surface 
Y = SOPRE(X,Z), then the computed surface is replaced by 
SOPRE(X,Z) at such points. The procedure is repeated with more 
line relaxations until the computed and assigned distributions 
agree. The surface modifications are determined from a first 
order partial differential equation that has been discussed 
earlier. Parameters controlling the scheme are discussed in 
the glossary. Assigning unrealistically low velocities near 
the leading and trailing edges serves to drive the computed 
surface onto the prescribed surface, which provides trailing 
edge closure and leaves the nose unaffected. 

The design distribution is defined by two or more section 
Mach number distributions. Linear interpolation is used to 
specify the values elsewhere. The section speed distributions 
are assigned over the computational domain. For a fixed cross 
section Z the lower trailing edge appears on the extreme left, 
the leading edge appears near the center and has the largest 
values of SO, and the upper trailing edge appears on the extreme 
right. 

The section distributions are defined by specifying input 
speeds Ql, Q2, Q3 and Q4 at fractions PCQl, PCQ 2, PCQ3 and 
PCQ4 =1 of the local chord and by interpolating in between 
with an exponential spline. The spline has free parameters 
that can be adjusted to prevent unwanted oscillations that would 
occur if cubic splines were used. In addition, a weighting 
parameter gives sagging curvature to the distribution so that 
two-dimensional shockless distributions can be simulated. 

The value Ql at the nose should be set so that the result- 
ing distribution lies below the initial analysis distribution 
along the lower surface and leading edge regions. Similarly, 
the value of Q4, the speed at the trailing edge, should be lower 
than the corresponding value computed in the analysis run. The 
two intermediate values, 02 and Q3, define the size of the 
supersonic zone and the section lift. The prescribed wing 
surface can be thinned out near the supersonic zone by removing 
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material smoothly to produce a slight groove. The depth and 
extent of the groove are determined by the three parameters 
DSURF, PCSl and PCS2. 

When the optimization routine is used a sequence of 
calculations is performed in the inverse mode to determine 
the gradient of the wave drag with respect to the assigned Mach 
numbers that define the design distribution. After the gradient 
is obtained, a line search of five steps is performed to mini- 
mize the drag. This procedure can be adjusted by changing the 
parameters that appear in subroutines OPT and VAIOA. 

The graphic output is produced in subroutines THREED and 
DRAGC and at the end of the main routine FL22INV. These programs 
have been written for the CDC 6600 at the Courant Institute 
of Mathematical Sciences and should be replaced by the plotting 
routines used at local installations or by dummy subroutines 
with the same names so that runs can be made without graphics . 

Output appears in both printed and graphical form. All 
the input data is immediately printed as output so that it is 
easy to check the accuracy of the input. 

At the beginning the coordinates defining the first span 
station are printed. If all the sections are similar only the 
coordinates of the leading edge, the chord and the twist are 
printed at the other stations. If the sections are different 
then the corresponding input profiles will be printed. The 
program prints the coordinates of the unfolded sections produced 
by the square root transformation at the root and the tip. 

A two-dimensional chart of the plane Y = 0 is printed giving 
values of an indicator IV which shows the properties of points 
in this coordinate surface. IV = 2 specifies a point on the 
wing; IV = 1 specifies a point on the trailing vortex sheet; 

IV = 0 specifies a point on the singular line X = 0; IV = -1 
specifies a point adjacent to the wing or vortex sheet; and 
IV = -2 specifies a point beyond the wing or vortex sheet. 

The iteration history is printed next. The maximum correc- 
tion to the velocity potential and the maximum residual of the 
difference equations together with its i,j,k location are 
printed at every cycle. For an analysis run the lift coeffi- 
cient CL, the wave drag coefficient CDW, two relaxation factors 
PIO and P20, a convergence factor BETA, and the number of 
supersonic points are printed at every iteration. For a design 
run, in addition to the correction and residual, the average 
difference between the computed and desired speeds and the 
corresponding maximum difference together with its i,k location 
are printed. The iteration cycles terminate after a given 
number of iterations or after a convergence criterion has been 
satisfied. A chart is then printed of the wave drag at the 
grid points (X,Z) of the wing surface. Supersonic points are 
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indicated by drag numbers IDRA.G ^ 0 and subsonic points are 
indicated by -1. 

After this, results for each span section are displayed. 

The section lift, drag and moment coefficients are printed. 

For an analysis run the mapped wing surface and the Mach 
niamber distribution are displayed as a printer plot. For a 
design run, the prescribed surface and the computed surface 
are shown. The prescribed and computed Mach number distribu- 
tions along the chord are shown for each span section. There 
are also Calcomp plots. The upper surface pressure distribu- 
tion at each span section and the corresponding wing sections 
with markings that indicate the wave drag on shocks are plotted. 
In an analysis run the same plots occur for each mesh refinement. 

The final printed results are the characteristics of the 
wing as a whole. These include the coefficients of lift, form 
drag, friction drag, total drag, the ratios of lift to form 
drag and lift to total drag, the pitching, rolling and yawing 
moments , and the wave drag . 

For an analysis run, the program repeats the same sequence 
of calculations and output on successively refined meshes. 


GLOSSARIES 


The input files consist of sequences of pairs of cards. 

The first card of each pair gives the names of the parameters 
that appear on the data card that follows. Each data card 
contains up to eight parameters with 10 colximns provided for 
each. The first input file described is needed for both analysis 
and design. If the code is used for analysis we are concerned 
with Tape 5 only and Card Pair 3 below does not exist. The 
input parameters are given in the order of their appearance on 
the input file. The input data is given in floating point 
format. The integer parameters are converted to integers in 
the code . 


Glossary of Tape 5 Parameters 


Card Pair 1: 

NX The number of mesh intervals in the direction 

of the chord. NX = 0 causes termination of the 
computation. 

NY The number of mesh intervals in the direction 

normal to the chord and span. 
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NZ 

FPLOT 

XSCAL, PSCAL 
FCONT 


F SWEEP 

Card Pair 2: 
MIT 

COV 

PIO 


The number of mesh intervals in the span 
direction. 

Controls the plots. 

FPLOT = 0. produces a print plot but no 
Calcomp plot. 

FPLOT ^ 1. produces a print plot and a 
Calcomp plot. 

These control the scales of the Calcomp plots. 

XSCAL = 0. scales each section plot to 5. 

PSCAL = 0. scales the pressure plots to 1. per 
inch. 

Determines the manner of starting the program. 

FCONT = 0. begins the calculation at itera- 
tion zero. 

FCONT = 1. continues the calculation from a 
previous calculation. For a design run the 
flow data (velocity potential and circulation) 
from an analysis run are read in on Tape 7 and 
used for initialization. The iteration count 
starts from zero for a design run. 

An indicator which selects the subroutine YST'TEEP 
used to solve the finite difference equations 
for the reduced potential G. 

FSWEEP =0. selects a vectorized YSWEEP 
subroutine . 

FSWEEP = 1. selects an unvectorized YSWEEP 
subroutine . 


The maximum number of iteration cycles which 
will be computed. 

The desired accuracy. If the maximum correc- 
tion is less than COV the calculation terminates 
or proceeds to a finer mesh. 

The subsonic relaxation factor for the velocity 
potential. PIO lies between 1. and 2. and 
should be increased linearly toward 2. with mesh 
refinement. 
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P20 The supersonic relaxation factor for the 

velocity potential. 

P20 1. Recommended value 1. 

P30 The relaxation factor for the circulation. 

Recommended value 1 . 

BETA The damping factor which controls the amount 

of added 4>st • Recommended value between 0. 
and 0.25. 

FHALF Determines whether the mesh will be refined. 

FHALF = 0. terminates the computation after 
MIT iterations or after convergence. 

FHALF ^ 0. halves the mesh after MIT cycles 
have been run on the crude mesh. An additional 
Card Pair 2 is required for each mesh refine- 
ment. The value FHALF = 0. appears on the last 
mesh refinement card. 

NDES Gives the number of surface modifications to be 

calculated. 

NDES £ 0. calls for an analysis run. 

NDES > 0. makes a design calculation with NDES 
surface modifications. MIT cycles of line 
relaxation are performed after each surface 
modification. No mesh refinements are made 
after the NDES cycles are completed. MIT = 1. 
with NDES ^ 100.- is recommended. If NDES > 0. 
an additional Card Pair 3 is needed at this 
point. 

Card Pair 3 : 

TSTEP TSTEP times the mesh increment in the X direc- 

tion is the time step defining the free boundary 
iteration. The recommended value is 0.03. 

FOO The coefficient multiplying the first order time 

derivative in the free boundary equation. This 
term dominates for subsonic flow. Recommended 
value FOO = 0.16. 

FIO The coefficient of the second derivative with 

respect to time and the X coordinate in the 
free boundary equation. This term controls 
convergence for supersonic flow. Recommended 
value FIO = -1. 
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FOPT 


FOPT ^ 0. indicates a regular inverse run. 
FOPT > 0. invokes the optimization procedure. 


Card Pair 4 : 

FMACH 

YAW 

ALPHA 

CDO 


The free stream Mach number. 

The yaw angle of the wing in degrees. 

The angle of attack in degrees. 

The estimated drag due to skin friction. 
This can be read in and added to the drag 
calculated by the program to give the total 
drag. 


Card Pair 5 : 
ZSYM 


NC 


SWEEPl 

SWEEP2 

SWEEP3 

DIHEDl 

DIHED2 


DIKED 


Indicates whether an isolated wing or a wing on 
a wall is being considered. 

ZSYM = 0. specifies an isolated wing at a 
prescribed yaw angle; obsolescent. 

ZSYM = 1. specifies a swept wing on a wall. 

The number of span stations from the wing root 
to the tip at which the wing section is defined 
if ZSYM = 1. For ZSYM = 0. the span stations 
are distributed from the leading to the trail- 
ing tip. The wing sections are each defined 
on subsequent cards. 

Sweep of the singular line at the wing root 
if ZSYM =1. or at the leading tip if ZSYM = 0. 

Sweep of the singular line at the tip. SWEEPl 
and SWEEP2 are used as the end conditions for 
the spline fit for the X coordinates of the 
singular line. 

Sweep of the singular line in the far field. 

Dihedral angle of the singular line at the wing 
root if ZSYM = 1. or at the leading tip if 
ZSYM = 0. 

Dihedral angle of the singular line at the 
wing tip. DIHEDl and DIHED2 are used as the 
end conditions for the spline fit of the 
Y coordinates of the singular line. 

Dihedral angle of the singular line in the far 
field. 
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Card Pair 6: 
Z 

XLE, YLE 
CHORD 

THICK 

ALPHA 

FSEC 


Card Pair 7: 
YSYM 


NU 

NL 

Card Pair 8: 
TRAIL 


Span location of the section. 

X and Y coordinates of the leading edge. 

The local chord value by which the profile 
coordinates are scaled. 

Modifies the section thickness. The Y 
coordinates are multiplied by THICK. 

The angle through which a section is rotated to 
introduce twist. 

Indicates whether or not the geometry for a new 
profile is supplied. 

FSEC = 0. means the section is obtained by 
scaling the profile used at the previous span 
section according to the parameters CHORD, 
THICK, and ALPHA. No further cards are read 
for this span station and the next card is the 
title card for the next span station, if any. 

FSEC = 1. means the coordinates for a new 
profile are to be read from the data cards 
that follow. 


Indicates the type or profile. 

YSYM = 0. means the data supplied are for a 
cambered profile. Coordinates are given for 
the upper and lower surfaces, each ordered 
from nose to tail with the leading edge 
included in both surfaces . 

YSYM = 1. means the data supplied are for a 
symmetric profile. A table of coordinates is 
read in for the upper surface only. 

The number of upper surface coordinates. 

The number of lower surface coordinates. For 
YSYM = 1 . , NL = NU . 


The included angle at the trailing edge in 
degrees. If the profile is open then TRAIL is 
the difference between the upper and lower 
trailing edge angles. 
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SLOPT 


XSING, YSING 


Card Pair 9: 
X, Y 

Card Pair 10: 
X, Y 


Card Pairs 11 


The last card 


The slope of the mean camber line at the 
trailing edge. This is used to continue 
the coordinate surface, assumed to contain 
the vortex sheet, smoothly off the trailing 
edge. 

The coordinates of the singular point inside 
the nose about which the square root trans- 
formation is applied to generate parabolic 
coordinates. This point should be located as 
symmetrically as possible between the upper 
and lower surfaces at a distance from the nose 
roughly proportional to the leading edge 
radius. The coordinates of the mapped profile 
in the output will show if this point has 
been located correctly. The coordinates of 
the singular point are chosen relative to the 
profile coordinates supplied in the cards which 
follow. 

(Upper surface coordinates.) 

The coordinates of the upper surface. They 
appear, a pair on each card, from leading edge 
to trailing edge. The format is (2F10.6). 

(Lower surface coordinates.) 

The coordinates of the lower surface from lead- 
ing edge to trailing edge. The leading edge 
point of the upper surface is the same as the 
leading edge point of the lower surface. The 
trailing edge points are different if the 
profile has an open tail. 


12,...: (Geometry at other span stations.) 

These cafds are like Card Pair 6, which defines 
Z, XLE, YLE, CHORD, THICK, ALPHA and FSEC for 
each section. The number of such cards depends 
on the number of span stations, NC. If FSEC = 0. 
new coordinates X,Y are not needed to define the 
profile . 

pair: 

The card which terminates the run is a repeat 
of Card Pair 1 with all the data set equal to 
zero . 
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Glossary of Tape 9 Parameters 


Card Pair 1: 
NQSTA 

Card Pair 2 : 
ZQSTA 

PCQl 

PCQ2 

PCQ3 

Q1 

Q2 

Q3 

04 

Card Pair 3; 
PCSl 

PCS 2 
DSURF 


The number of span stations at which the 
design distribution is defined from wing root 
to tip. 


(Card Pairs 2 and 3 are repeated NQSTA times . ) 

The Z coordinate of the span section at which 
the design distribution is given. 

The location of the first specified value Q1 
of the speed, expressed as a fraction of the 
local chord. 

Location of the second specified value Q2. 

Location of the third specified value Q3 . 

(PCQ4 = 1 because the speed is always 
prescribed at the trailing edge.) 

The first prescribed Mach number near the 
leading edge used in spline fitting the design 
distribution at each section. 


The prescribed Mach number near the front of 
the supersonic zone. 

The prescribed Mach number near the rear of 
the supersonic zone. 

The prescribed Mach niimber at the trailing edge. 


(These parameters are used to define the groove 
for each span station.) 


Location of the left edge of the groove 
expressed as a fraction of the local chord. 

The groove is assumed to be on the upper surface. 


Location of the right edge of the groove 
expressed as a fraction of the local chord. 

The maximum depth of the groove expressed in 
units used in the computational domain. 
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LISTING OF A SAMPLE RUN 


FPLOT 

0 . 


FNX 

. 120OOE+O3 
XSCAL 

. 56000E+01 


FIT(NM) 

. lOOOOE+01 


FNY 

.lOOOOE+02 

PSCAL 

50000E+00 


CQVO(NM) 

. lOOOOE-06 


FNZ 

.12000E+02 

FCONT 
. lOOOOE+01 


PIO(NM) 

.17200E+01 


FSWEEP 

.lOOOOE+Oi 


P20(NM) 

. lOOOOE+01 


P30(NM) BETAO(NM) FHALF(NM) FDES(NM) 

.lOOOOE+01 .lOCOOE+00 0. .lOOOOE+03 


TSTEP FOO FIO Fll 

.30000E-01 .16000E+00 -.lOOOOE+01 0. 

FQPT 

0 . 


FMACH YA 

.83000E+00 0, 


AL COO 

.lOOOOE+01 0. 


ZSYM 

. lOOOOE+01 


FNC 

.60000E+OI 

DIHEDl 

0 . 


XL 

0 . 


SWEEPl 

.30000E+02 

DIHED2 

0 . 


YL 

0 . 


SWEEP2 
. 30000E+02 

DIHED 

0 . 


SWEEP 

. 30000E+02 


ZS (K) 

0 . 


CHORD THICK AL FSEC 

.67370E+00 .lOOOOE+01 0. 0. 


YSYM 

0 . 


FNU FNL 

.A7000E+02 .35000E+02 


TRL 

.1577E+01 


SLT 

-.lOOOE+00 


XSING 

.lOOOE-01 


YSING 
. 1A03E-01 



UPPER SURFACE 


LOWER SURFACE 


XP ( I) 

.28010E-02 
.46130E-02 
.19370E-01 
.26231E-01 
.35059E-01 
.45<tl5E-01 
..57105E-01 
.69832E-01 
.83659E-01 
.98782E-01 
. 11506E+00 
.13213E+00 
.1^999E+00 
.16893E+00 
.1089faE+OO 
.20982E+00 
. 231-^66 + 00 
.25A00E+00 
.27723E+00 
.30076E+00 
.32450E+00 
. 34053E+OO 
.37262E+00 
. 396A5E+00 
.41996E+00 
. 94316E+00 
. A6569E+00 
. A8715E+00 
.50725E+00 
. 52567E+0C 
.!?4206E + 00 
. 55629E+00 
.56854E+00 
.57908E+00 
. 58821E+00 
. 69637E+00 
.60433E+00 
.61266E+00 
.67205E+00 
.73993E+00 
.80171E+00 
.86178E+00 
.91109E+00 
.95185E+00 
.97857E+00 
,99^11E+00 
.99673E+00 


YP( I) 

.22664E-01 
. 35103E-01 
.-«f6090E-01 
.^9082E“01 
.52122E-01 
.54992E-01 
.57681E-01 
. 60167E-01 
.62^97E-01 
.6A716E-01 
.66806E-01 
.68731E-01 
.70501E-01 
.721A8E-01 
.73667E-01 
.75033E-01 
.76243E-01 
.77299E-01 
. 78188E-01 
.78895E-01 
.79252E-01 
.79521E-01 
.79618E-01 
.79550E-01 
.79320E-01 
. 78932E-01 
.78400E-01 
,777^7E-01 
.76996E-01 
. 7618AE-01 
.753^6E-01 
.7A526E-01 
.73734E-01 
.72982E-01 
.72270E-01 
.71577E-01 
.70852E-01 
.70035E-01 
. 6267iE-Ql 
.51620E-01 
.39381E-01 
.26501E-01 
. 16163E-01 
.67060E~02 
.A7780E-U2 
.29260E-02 
.26110E-02 


XP( I) 

-.28010E-02 

0 . 

.80950E-02 
.19227E-01 
.33397E-01 
.49561E-01 
.67610E-01 
.87596E-01 
. 10960E+00 
.13392E+00 
.15988E+00 
.187A2E+00 
. 21663E+00 
.2A729E+00 
.27937E+00 
. 312fc^E+00 
.3^718E+00 
.38259E+00 
. A1913E+00 
.456366+00 
.494726+00 
. 53380E + 00 
. 57433E+00 
. 61612E+00 
.66026E+00 
. 70633E+00 
.755C4E+00 
.80455E+00 
.85410E+00 
. 89932E+00 
.93856E+00 
.967896+00 
.98810E+00 
.99652E+00 
.99733E+00 


YP(I) 

.22664E-01 

0 . 

-.98270E-02 

-.172556-01 

-.23432E-01 

-.285726-01 

-.328916-01 

-.365896-01 

-.39868E-01 

-.42776E-01 

-.453396-01 

-.47551E-01 

-.49431E-01 

-.50965E-01 

-.52158E-01 

-.528006-01 

-.531456-01 

-.53076E-01 

-.525376-01 

-.51470E-01 

-.498086-01 

-.47365E-01 

-.439576-01 

-.39376E-01 

-.33304E-01 

-.25802E-01 

-.172816-01 

-.926306-02 

-.35060E-02 

-.109206-02 

-.13510E-C2 

-.273106-02 

-.415706-02 

-.49110E-02 

-.50020E-02 
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SECTION DEFINITION AT Z » 0.00000 

XLE YL£ CHORD THICKNESS RATIO 

0.0000 0.0000 .6737 1.0000 

ZS(K) XL YL 

.20000E+00 .11500E+00 0. 

CHORD THICK AL FSEC 

.61470E+00 .lOOOOE+01 0. 0. 

SECTION DEFINITION AT Z » .20000 

XLE YLE CHORD THICKNESS RATIO 

.1150 0.0000 .6147 1.0000 

ZS(K) XL YL 

.40000E+00 .23000E+00 0. 


CHORD THICK AL FSEC 

.55580E+00 .lOOOOE+01 0. 0. 

SECTION DEFINITION AT Z • .40000 

XLE YLE CHORD THICKNESS RATIO 

.2300 O.OOOC .5558 1.0000 

ZS(K) XL YL 

.60000E+00 .34500E+00 0. 

CHORD THICK AL FSEC 

.49680E+00 .lOOOOE+01 0. 0. 

SECTION DEFINITION AT Z » .60000 

XLE YLE CHORD THICKNESS RATIO 

.3450 0.0000 .4968 1.0000 

ZS(K) XL YL 

.80000E+00 .46000E+00 0. 


CHORD THICK AL FSEC 

.43790E+00 .lOOOOE+01 0. 0. 

SECTION DEFINITION AT Z ■ .60000 

XLE YLE CHORD THICKNESS RATIO 

.4600 0.0000 .4379 1.0000 

ZS<K) XL YL 

.lOOOOE+01 .57500E+00 0. 


CHORD THICK AL FSEC 

.378906+00 .lOOOOE+Ol 0. 0. 

SECTION DEFINITION AT Z « 1.00000 

XLE YLE CHORD THICKNESS RATIO 

.5750 0.0000 .3789 1.0000 


ALPHA 

0.0000 


ALPHA 

0.0000 


ALPHA 

0.0000 


ALPHA 

0.0000 


ALPHA 

O.OGOO 


ALPHA 

0.0000 
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vj> u) ro !-■ 


PARAMETERS TO DEFINE THE ASSIGNED DESIGN MACH NUMBER DISTRIBUTION 


z 

PCMl 

PCM2 

PCM3 

Ml 

M2 

M3 

M4 

0.000 

.065 

.220 

.830 

.420 

.650 

1.240 

.670 

.250 

.055 

.220 

.800 

.440 

1.140 

1.190 

.660 

.500 

.055 

.220 

.780 

.440 

1.270 

1.210 

.670 

.875 

.065 

.220 

.780 

.480 

1.370 

1.120 

.660 

1.000 

.065 

.200 

.820 

.500 

1.345 

1.000 

.700 


PCXl 

PCX2 

OSURF 

. 500 

.900 

.002 

. 400 

.900 

.002 

. 400 

.900 

.002 

. 150 

.900 

0.000 

. 150 

.900 

0.000 
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INDICATION IV ( I ^ K ) OF WING AND VORTEX SHEET IN PLANE Y «0 
( ( IV ( I # K ) iK » Kl > K 2 )> I « 2 > NX ) 
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2 2 2 2 2 
2 2 2 2 2 
2 2 2 2 1 
2 2 2 2 1 
2 2 2 1 1 
2 2 2 1 1 
2 2 111 
2 1111 
2 1111 
11111 
11111 
11111 
11111 
11111 
11111 
11111 
11111 
11111 
11111 


2 1 
1 1 
1 1 
1 1 
1 1 
1 1 
1 1 
1 1 
1 1 
1 1 
1 1 
1 1 
1 1 
1 1 
1 1 
1 1 
1 1 
1 1 
1 1 


1 1 - 1-2 -2 

1 1 - 1-2 -2 

1 1 - 1-2 -2 

1 1 - 1-2 -2 

1 1 - 1-2 -2 

1 1 - 1-2 -2 

1 1 - 1-2 -2 

1 1 - 1-2 -2 

1 1 - 1-2 -2 

1 1 - 1-2 -2 

1 1 - 1-2 -2 

1 1 - 1-2 -2 

1 1 - 1-2 -2 

1 1 - 1-2 -2 

1 1 - 1-2 -2 

1 1 - 1-2 -2 

1 1 - 1-2 -2 

1 1 - 1-2 -2 

1 1 - 1-2 -2 
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NORMAL CELL DISTRIBUTION IN SQUARE ROOT PLANE 


Y 

.56695 
.44721 
.35909 
. 28868 
.22917 
.17678 
.12910 
.08452 
.04181 
0.00000 

SCALE FACTOR POWER LAW 
.50000 .50000 


SPANWI 


0 . 


1 . 

1 . 

1 . 

1 . 


TIP 


SE CELL DISTRIBUTION AND 


z 

X SING 

00000 

. 00862 

12500 

.08018 

2500C 

.15138 

37500 

.22282 

50000 

.29424 

62500 

.36565 

75000 

.43709 

87500 

. 50829 

00000 

.57985 

12677 

. 65304 

26516 

.73294 

43301 

.82985 


YZ 

. 00000 
.00333 
.00271 
.00233 
.00261 
.00233 
.00271 
.00333 
.00000 
0.00000 
0.00000 
0.00000 

LOCATION POWER LAW 
57143 .50000 


IN6ULAR LINE 


Y SING 

XZ 

-.00582 

. 57735 

-.00556 

. 56934 

-.00516 

. 57083 

-.00486 

.57173 

-.00454 

.57106 

-. 00423 

. 57173 

-.00393 

.57083 

-.00353 

. 56934 

-.00327 

.57735 

-.00327 

. 57735 

-.00327 

. 57735 

-.00327 

. 57735 


XZZ 

YZZ 

-. 10644 

.04426 

-.02173 

.00902 

.01938 

-.00807 

-.00491 

.00206 

-.00000 

. 00000 

. 00491 

-.00206 

-.01938 

.00807 

.02173 

-.00902 

.10644 

-.04426 

0.00000 

0.00000 

0.00000 

0.00000 

0. 00000 

0.00000 
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ITERATIVE SOLUTION 

STRIP WIDTH FOR HORIZONTAL LINE RELAXATION 

1.00000 

NX NY NZ 

128 10 12 


COMPUTING TIME A7.497 SECONDS 


MACH NO 

YAW 


ANG 

OF 

ATTACK 


.83000 

0.00000 


1.00000 


ITERATION 

RESIDUAL 

I 

J 

K 

AVERAGE a 

DRAG 

1 

.32652E-05 

128 

10 

6 

0. 

.0028 

2 

.22339E-04 

128 

10 

5 

.24099E+00 

.0026 

3 

-.33525E-04 

110 

11 

4 

.23144E+00 

.0028 

A 

-.66452E-04 

ill 

11 

3 

.22337E+C0 

.0027 

5 

-.65743E-04 

111 

11 

3 

.21972E+00 

.0026 

6 

-.61413E-04 

111 

11 

3 

.21434E+00 

.0025 

7 

-.55446E-04 

111 

11 

3 

.20744E+00 

.0024 

8 

-.49685E-04 

111 

11 

3 

.20334E+00 

.0023 

9 

-.42556E-04 

111 

11 

3 

.19688E+00 

.0022 

10 

-. 39147E-04 

111 

11 

3 

.19178E+00 

.0022 

11 

.39657E-04 

106 

11 

3 

. 18661E+00 

.0021 

12 

.40238E-04 

106 

11 

3 

.18115E+00 

.0020 

13 

.40815E-04 

106 

11 

3 

. 17402E+00 

.0020 

14 

. 41442E-04 

106 

11 

3 

.16688E+00 

.0019 

15 

.41479E-04 

106 

11 

3 

. 16052E+00 

.0019 

16 

.42238E-04 

118 

11 

3 

.15602E+00 

.0018 

17 

.44582E-04 

118 

11 

3 

. 15039E + 00 

.0018 

18 

.46007E-04 

118 

11 

3 

.14704E+00 

.0017 

19 

.47812E-04 

118 

11 

3 

.14158E+00 

.0017 

20 

. 46790E-04 

118 

11 

3 

.13557E+00 

.0016 

21 

.49760E-04 

118 

11 

3 

. 13093E+00 

.0016 

22 

.50483E-U4 

118 

11 

3 

.12642E+00 

. 0016 

23 

. 50944E-04 

118 

11 

3 

.12177E+00 

.0015 

24 

. 51533E-04 

118 

11 

3 

. 11737E + 00 

.0015 

25 

. 5160CE-04 

118 

11 

3 

. 11316E + 00 

.0015 

26 

.51489E-04 

118 

11 

3 

.10941E+00 

.0015 

27 

. 51818E-04 

118 

11 

3 

.10650E+00 

.0015 

28 

.51418E-04 

118 

11 

3 

. 10299E+00 

.0014 

29 

.51416E-04 

118 

11 

3 

.99338E-01 

.0014 

30 

.51081E-04 

118 

11 

3 

.96995E-01 

.0014 

31 

. 50849E-04 

118 

11 

3 

.93886E-01 

.0014 

32 

.50622E-04 

118 

11 

3 

.91089E-01 

.0014 

33 

.50489E-04 

118 

11 

3 

. 88611E-01 

.0014 

34 

. 50499E-04 

118 

11 

3 

. 86250E-01 

.0014 

35 

.50617E-04 

118 

11 

3 

.82413E-01 

.0014 

36 

. 50744E-04 

118 

11 

3 

.79798E-01 

.0014 

37 

.50766E-04 

118 

11 

3 

. 77769E-01 

.0014 

38 

. 50833E-04 

118 

11 

3 

. 76239E-01 

.0014 

39 

.50722E-04 

118 

11 

3 

.74098E-01 

.0014 

40 

. 50504E-04 

118 

11 

3 

.71992E-01 

.0014 

41 

.50140E-04 

118 

11 

3 

.70111E-01 

.0014 

42 

. 49725E-04 

118 

11 

3 

. 68546E-01 

.0014 

43 

.49258E-04 

118 

11 

3 

.66656E-01 

.0014 

44 

.48772E-04 

118 

11 

3 

.646738-01 

.0014 

45 

.48218E-04 

118 

11 

3 

. 63204E-01 

.0014 

46 

.47602E-04 

118 

11 

3 

.62074E-01 

.0014 



47 

.46976E-04 

118 

11 

3 

. 60766E-01 

.0014 

48 

.46369E-04 

118 

11 

3 

.59644E-01 

.0014 

49 

.45817E-04 

118 

11 

3 

.58050E-01 

.0014 

50 

.45222E-04 

118 

11 

3 

.56752E-01 

.0014 

51 

.44648E-04 

118 

11 

3 

.55776E-01 

.0014 

52 

.44038E-04 

118 

11 

3 

.54843E-01 

.0014 

53 

. 43466E“04 

118 

11 

3 

.53823E-01 

.0014 

54 

.42837E-04 

118 

11 

3 

. 52818E-01 

.0013 

55 

.42229E-04 

118 

11 

3 

.51770E-01 

.0013 

56 

.41618E-04 

118 

11 

3 

. 50704E-01 

.0013 

57 

.41004E-04 

118 

11 

3 

.49794E-01 

.0013 

58 

.40368E-04 

118 

11 

3 

.48961E-01 

.0013 

59 

.39772E-04 

118 

11 

3 

.477938-01 

.0013 

60 

.39209E-04 

118 

11 

3 

.46724E-01 

.0013 

61 

.38655E-04 

118 

11 

3 

.45672E-01 

.0013 

62 

.38039E-04 

118 

11 

3 

.45067E-01 

.0012 

63 

.37509E-04 

118 

11 

3 

.44301E-01 

.0012 

64 

.36847E-04 

118 

11 

3 

. 43454E-01 

.0012 

65 

.36252E-04 

118 

11 

3 

.42829E-01 

. 0012 

66 

.35627E-04 

118 

11 

3 

.42009E-01 

.0012 

67 

.35008E-04 

118 

11 

3 

. 41134E-01 

.0012 

68 

. 34363E-04 

118 

11 

3 

. 40542E-01 

.0012 

69 

.33749E-04 

118 

11 

3 

. 39980E-01 

.0011 

70 

.33131E-04 

118 

11 

3 

. 39322E-01 

.0011 

71 

.32486E-04 

118 

11 

3 

. 38786E-01 

.0011 

72 

.31794E-04 

118 

11 

3 

.37994E-01 

.0011 

73 

.31138E-04 

118 

11 

3 

. 37112E-01 

.0011 

74 

.30488E-04 

118 

11 

3 

. 36548E-01 

.0011 

75 

.29836E-04 

118 

11 

3 

. 35976E-01 

.0011 

76 

.29193E-04 

118 

11 

3 

.35314E-01 

.0010 

77 

.28521E-04 

118 

11 

3 

.34710E-01 

.0010 

78 

.27883E-04 

118 

11 

3 

.34138E-01 

.0010 

79 

.27184E-04 

118 

11 

3 

.33668E-01 

.0010 

80 

.26540E-04 

118 

11 

3 

. 33157E-C1 

.0010 

81 

.25852E-04 

118 

11 

3 

.325438-01 

.0010 

82 

.25194E-04 

118 

11 

3 

.32108E-01 

.0010 

83 

.24503E-04 

118 

11 

3 

. 31485E-01 

.0009 

84 

.23841E-04 

118 

11 

3 

.31071E-C1 

.0009 

85 

. 23131E-04 

118 

11 

3 

.30416E-01 

.0009 

86 

.22452E-04 

118 

11 

3 

.29559E-01 

.0009 

87 

.21766E-04 

118 

11 

3 

.29395E-01 

.0009 

38 

.21059E-04 

118 

11 

3 

. 29200E-01 

.0009 

89 

.20368E-04 

118 

11 

3 

.28790E-01 

.0009 

90 

.19682E-04 

118 

11 

3 

.28600E-01 

.0009 

91 

.18975E-04 

118 

11 

3 

.28460E-01 

.0009 

92 

.18292E-04 

118 

11 

3 

.28356E-01 

.0008 

93 

.17604E-04 

118 

11 

3 

.28030E-01 

.0008 

94 

.16921E-04 

118 

11 

3 

.27879E-01 

.0008 

95 

.16242E-04 

118 

11 

3 

.27734E-01 

.0008 

96 

. 15568E-04 

118 

11 

3 

.27594E-01 

.0006 

97 

.14893E-04 

118 

11 

3 

.27635E-01 

.0008 

98 

. 14235E-Q4 

118 

11 

3 

.27428E-01 

.0008 

99 

.13593E-04 

118 

11 

3 

.27197E-01 

.0008 

100 

.12961E-04 

118 

11 

3 

.26965E-01 

.0008 

MAX RESIOAL 

1 MAX RESIOAL 2 


WORK 

REDUCTN/CrCLE 


.3265E-05 ,12y6E-0^ 99.0000 1.01^0 

COMPUTING TIME 672.962 SECONDS 
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WAVE DRAG « .00078 

PRINTOUT OF IDRAGd/K) 
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0 
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0 

0 

0 

0 

0 

6 

-1 

0 

0 

0 

0 

0 

0 

4 

23 

-1 

0 

0 

0 

0 

0 

0 

25 

-1 

-1 

0 

0 

0 

0 

0 

10 

22 

-1 

-1 

0 

0 

0 

0 

13 

57 

“1 

-1 

-1 

0 

0 

0 

0 

75 

19 

-1 

-1 

“1 

0 

0 

0 

15 

25 

-1 

-1 

-1 

-1 

0 

0 

0 

66 

-1 

-1 

-1 

-1 

-1 

0 

0 

23 

21 

-1 

-1 

-1 

-1 


0 

0 

59 

-1 

-1 

-1 

-1 

-i 


0 

66 

20 

-1 

-1 

-1 

“1 



0 

149 

-1 

-1 

-1 

-1 

-1 



19 

17 

-1 

-1 

-1 

-1 




90 

-1 

-1 

-1 

-1 

-1 




47 

-1 

-1 

-1 

-1 





-1 

“1 

-1 
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-1 
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-1 

-1 

-1 







-1 
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-1 -1 
-1 
-1 
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WING CHARACTERISTICS 


MACH NO 
.630C0 

CL 

.40535 


CD WAVE 
.00078 

CM PITCH 
-.39938 


YAW 

ANG OF ATTACK 


0.00000 

l.OOCOO 


CD FORM 

CD FRICTION 

CO 

.01729 

O.COCOO 

.01729 

L/D FORM 

L/D 


23.439b6 

23.43966 



CM ROLL 

CM YAW 

AWIN 

.37631 

.00116 

.52764 
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LISTING OF THE CODE 


PROGRAM F L22 1 NV ( INPUT»512> OUTPUT-512# TAPES- INPUTS T AP E 6-OUT PUT # T AP E 
114-0# TAPE 15-0# TAPE 16-0# TAPE 1-0# TAPE?- 51 2#T APE b«5 12# TAP Ell- 0# TAPE 9- 
264#TAPE10-512) 

C MAIN ROUTINE WHICH CONTROLS THE COMPUTATIONAL PROCEDURE. 

C G IS REDUCED VELOCITY POTENTIAL 

COMMON G(129# 12#15)#S0(129#15 )#E0(131)#Z0( 131)# IV(129#15)# ITEK15) 
1# ITE2( 15)# AO ( 129)# Al( 129)# A2( 129)# A3 (129) #80(12 )# 81(12 )#B2 (12)# B3( 
2i2)#Z(15)#Cl(15)#C2(15)#C3(15)#XC(15)#XZ(15)#XZZ(15)#YC(15)#Y2(15) 
3#YZZ( 15)#NX#NY#NZ#KT£1#KTE2# I5YM#KSYM, SCAL# SCALZ# YAW# CYAW# SYAW# ALP 
4HA#CA#SA#FMACH#N1#N2#N3#I0#NDES#TSTEP#£PS1#UPRE ( 129# 15 ) # SOPR E ( 129# 
515)#NQSTA# ZQSTA( 15)#PCQ1(15)#PCQ2(15)#PCQ3(15)#QQ1( 15)#QQ2(15)#QQ3 
6(15)#QQ4(15)# PCSK 15)# PCS2(15 )# DSURF( 15) #RDQ#RDSO# FOO# FOl# FIO# Fll# 
7NDQ# IQ#KQ# AWING# VOL DRG# I DRGP LT ( 129# 15 ) # S EC DRG ( 1 5 ) 

COMMON /FLO/ STRIP#P1#P2#P3#BETA#FR#IR#JR#KR#DG#IG#JG#KG#NS#FSWEEP 
DIMENSION XS( 200#1D# YS(200#1D# ZS(ll)# XLE(ll)# YLE(ll)# SLOPT( 
111)# TRAIL(ll)# NP(1D# EKID# E2(ll)# E3(1D# E4(1D# E5(1D# XP 
2(241)# YP(24D# 01(241)# D2(24D# 03(241)# X(129)# Y(129)# SV(129) 
3# SM(129)# CP(129)# CH0RD(15)# SCU15)# SCD(15)# SCM(15)# FIT(3)# 
4CDV0(3)# P10(3)# P20(3)# P30(3)# BETA0(3)# STRIPCO)# FHALF(3)# RE 
5S(50D# CDUNT(50D# UC(129)# VC(129)# WC(129)# FDES(3)# CLU(ll)# C 
6LPRE(15)# ALF0(15) 

DIMENSION DESC(IO)# LABEL(IO)# NPARAM(30)# TITLE(IO) 

COMMON /DIM/ NX1#NY1#NZ1#FDIM 

ND-200 

NE-129 

IREAD-5 

IWRIT-6 

KPLOT-0 

IPLOT-1 

ISTQP-2 

JO-0 

NFl-1 

REWIND 7 

RAD-57,2957795130823 
10 WRITE ( IWRIT# 670) 

WRITE (IWRIT#39C) 

READ (IREAD#660) TITLE 
WRITE (IWRIT#700) TITLt 
READ (IREAD#660) DESC 

READ (IREAD#650) FNX# FNY# FNZ# FP LOT# XSCAL# PSCAL# FCONT# FSWEEP 

WRITE (IWRIT# 780) FNX# FN Y# FNZ# F P LOT# X SC AL# P SC AL # FCONT # F S WE EP 

NX-FNX 

NY-FNY 

NZ-FNZ 

IF (NX.LT.l) GO TO 360 
KPLOT-ABS( FPLOT) 
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RbAD (IREAD>660) DESC 
WRITE (IWRIT^790) 

NH-0 

20 NM-NM+1 

READ (I READ# 6 50) F I T ( NN ) # C OVO ( NM ) # P 10 ( NM ) # P20 ( NM ) # P 30 ( NM ) # DE TAO ( NM 
1)#FHALF(NM)#FDES{NM) 

STRIP0(MM)-1.0 

WRITE { I WRIT# 690) F IT ( NM ) # COVO( NM ) # P 10 ( NM ) # P20 ( NM ) # P30 (NM ) # BET A0( N 
iri)#FhALF(NM)#FDES(NM) 

IF (FHALFINM) .NE.C, .ANO.NM.lt. 3) GO TO 20 
IF ( FDESd ) .LE.O. ) GO TO 30 
READ (IREAD#660) DESC 

READ (IREA0#650) TS T£ P# FOO# FIO# F 1 1# F0P7 
WRITE (IWRIT#AOO) DESC 

WRITE (IWRIT#A10) TS TE P # F 00# F 10# Fll# FDP T 
30 CONTINUE 
NMESH-NM 
FHALF(3)«0. 

READ (IREAD#660) DESC 

READ (IREAD#650) FM AC H# Y A# A L# CDO# C LOP T# RC L# SR EF 
WRITE (IwRIT#800) FM AC H# Y A# AL#C DO# SRE F 
YAW-Y A/RAD 
ALPHA-AL/RAO 

CALL GEQM {ND#NC#NP#2S#XS#YS#XLE#YLE#SL0PT#TRAIL#XP#YP#SWEEP1#SWEE 
1P2#SWEEP#DIHE01#DIHED2#DIHED# XTEO#CHDRDO#ZTIP# ISYM0#KSYM#CL0) 
ISYM-ISYMO 

IF ( ALPHA. NE.O. ) ISYM-0 
IF (KSYM.NE.O) YAW-O. 

CYAWCQS ( YAW) 

SYAW*SIN ( YAW) 

CA»CYAW>^COS (ALPHA ) 

SA*CYAW^SIN( ALPHA) 

IF ( FDES ( 1 ) .GT.O. ) CALL READQS ( NQ S T A# ZQST A# PC Q 1 # PC Q2# PC Q3 # QQ 1 # QQ 2 
1# QQ3#UQA#PCS1#PCS2,DSURF#FMACH) 

IF (FCONT.lt. 1.) GO TO 50 

READ (7) NX#NY#NZ#NM#K1#K2#N1T 

MX-NX+1 

MY-NY+2 

MZ»NZ+3 

IF ( FDESd) .GT.O. ) NM»1 
IF (FDES(NM).GT.O.) MT-0 
DO 40 K»1#MZ 

READ (7) ( (G( I# J#K)# I-1#MX)# J-1#MY) 

40 CONTINUE 

READ (7) (E0(K)#K»K1#K2) 

REWIND 7 
50 CONTINUE 

FDIM-FHALF (1 ) 

NX1-NX+40-20+FDIM 

NYl»NY+2-FDIM 

NZl»NZ+2-l*FDIM 

CALL COORD (NX#NY#NZ#KSYM# XTEO# ZTIP# XMAX# ZM AX# SY# SC AL # SC AL Z# AX# AY# 
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lAZfAC^Ali A2^ A3f BC>Bl>d£^B3^Z^Cl/C2>C3) 

CALL SINGL ( N C, NZ # KS YM> KT E 1/ K TE 2^ CHOR DC # S W b EP 1# S kJE E P2 # S WE E 0 1 HED 1 
1#DIHED2#DIHED#2S>XLE>YLE>XC>XZ/XZZ, YC>YZ»YZZ^ Z>C1^C2#C3# fcl#E2#E3>E 
24>E5»IND>CL0#CLPRE) 

CALL SURF (ND#NE>NC,NX>NZ> I S YM> K SYM> KTE 1# K TE 2> SC A L# YAW> AO» Z»ZS>XC# 
lYCfSLOPTf TRAILS XS^YS^NP^ ITEl^ 1TE2» i V# SO# ZO^ XP# YP # D1 j D2» D3# X/ Y ^ IN D, 
2XZ#YZ#A1#C1) 

IF (IND.EQ.O) GO TO 370 

IF (FDES(l) .GT.O. ) CALL SETQS ( N E# N X# QP RE# SO# SOP R £# I T E 1# I T E2 # KT El # 
1KTE2#Z#ZQSTA# A0#PCQ1#PCQ2#PCU3#UC#VC#QC1#Q02#0Q3#Q04#PCS1#PCS2#DSU 
2RF#NQSTA) 

IF (FCONT.GE. 1. ) GO TO 60 

NM»1 

NIT-0 

CALL ESTIM (ALFO) 

60 WRITE (IWRIT#670) 

FCQNT-0. 

MIT»FIT(NM)+NIT 

KRES-2 

JRES=0 

NRES-0 

COV-CQVO(NM ) 

STRIP-STRIPO(NM) 

BETA«BETAO(NM ) 

MX-NX+1 

MY-NY+2 

MZ-NZ+3 

KY-NY+1 

Kl-2 

K2-NZ 

IF (KSYM.EQ.O) GO TO 70 
Kl-3 
K2-NZ+2 
70 LZ-NZ/2+1 

IF (KSYM.NE.O) LZ-3 
WRITE (IWRIT#A20) 

DO 80 I=2#NX 

80 WRITE (IWRIT#720) ( I V ( I # K ) # K -K 1 # K2 ) 

WRITE (IWRIT#670) 

WRITE (IWRIT#^30) 

DO 90 I-2#NX 

90 WRITE (IWRIT#680) AO ( I ) # SO ( I # L Z ) # SO ( I # K TE2 ) 

WRITE (IWRIT#AAO) 

WRITE (IWRIT#680) XNAX#AX 
WRITE (1WRIT#670) 

WRITE (IWRIT#450) 

DO 100 J-2#KY 

100 WRITE (IWRIT#680) BO(J) 

WRITE {IWRIT#A60) 

WRITE (IWRIT#680) SY#AY 
WRITE (IWRIT#670) 

WRITE (IWRIT#470) 
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DQ 110 K-K1#K2 

110 WRITE (IWRITfbSO) Z ( K ) ^ XC ( K ) / YC ( K ) « XZ (K ) > YZ (K ) > X ZZ ( K ) > YZZ ( K ) 
WRITE (IWRIT«460) 

WRITE (IWRIT>660) ZMAX/AZ 
WRITE (IWRIT»670) 

WRITE (IWRIT^490) 

WRITE (IWRIT>680) STRIP 
WRITE (IWRIT^SOO) 

WRITE (IWRIT>710) NX>NY^NZ 
CALL SECOND (T) 

WRITE (IWRIT>770) T 
WRITE (IWRIT#510) 

WRITE (IWRIT#680) FMACH>YA,AL 

IF (FOES(NH).LE.O..ANO.CLOPT.LE.O.) WRITE (IWRIT>5<rO) 

IF ( FOES(NM) .GT.O. ) WRITE <IWRIT#530) 

IF (CL0PT.6T.0. ) WRITE (IWRIT>520) 

KDES-0 

NDES-FDES (NM) 

LX-NX/2+1 

CL-0. 

DO 120 K-K1/K2 
Il-ITEKK) 

X(I1)-XC(K) + . 5*SCAL*{ A0(I1)*A0( I1)-S0(I1,K )*S0( I1#K ) ) 

X(LX) «XC (K)+. 5*SCAL*(A0(LX)*A0{ LX)-SO(LX,K )*SO(LX,K) ) 
CH0RD(K)-X(I1)-X(LX) 

120 CONTINUE 

KZDUM-KTE2-1 

S-0. 

DO 130 K-KTE1>KZDUM 
DZ0».5*(Z{K+1)-Z(K) ) 

130 S -S+DZO* ( CHORD (K+1 ) +CHORD ( K) ) 

AWIN6-S 

140 KDES-KDES+1 

IF (NDES.GT.O) NIT-0 
150 NIT-NIT+1 
Pl-PIO(NM) 

P2»P20(NM) 

P3-P30(NH) 

IF (FOPT.LT.l.) GO TO 160 
CALL OPT (QQ1>QQ2>QQ3>0Q4) 

GO TO 250 
160 CONTINUE 

CALL MIXFLO 
FCL-0. 

KCL-0 

IREFIN-0 

VOLDRG-0. 

CALL ORAGC (0#0.) 

IF (NDES.GT.O) 60 TO 180 
DO 170 K-3#MZ 
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IF (K.LT.KTE1.0R.K.GT.KTE2) GO TO 170 
CALL VELQ { K* K, S V# SM» C P # X# Y, UC > VC> WC ) 

Il-ITEKK) 

I2«ITE2(K) 

CH0RD(K)«X(I1)-X(LX) 

CALL FORCF ( 1 1 2^ X> Y/C P/ AL ^CHORD ( K > / XC ( K > / SCL ( K >/SCD(K)^5CM(K)) 
170 CONTINUE 

CALL TOT FOR ( KTE 1# KTE2# CHORD# SCL^SC D, SCM# 2/ XC #C L > CD1> CMP > CMR» CMY» A 
IWING) 

180 CONTINUE 
J0»0 

IF (NDES.LE.O) GO TO 190 
IF (NDQ.GT.O) RDQ»RDQ/FLOAT(NDQ ) 

IF (CLOPT.LE.O. ) WRITE (IWRIT#7<*0) KDES# D6# IG# J G# KG# FR, 1 R# JR, KR# RD 
1Q#RDS0#IQ#KQ# VOLDRG#CL 

190 IF (NDES.LE.O.AND.CLOPT.LE.O. ) WRITE (IWRIT#730) N I T# DG# I G# J G# KG# F 
1R#IR# JR#KR#CL#V0LDRG#P1#P2#BETA#NS 
IF (CLQPT.GT. 0. ) WRITE (IWRIT#750) N IT# DG# I G# JG # KG# FR# I R# J R# K R# FC L 
1#KCL#RCL#NS 
JRES» JRES+1 

IF ( JRES.EQ.KRES) JRES=1 
IF (JRES.NE.l) GO TO 200 
NRES»NRES+1 
C0UNT(NRES)»NIT-1 

IF (NDES.6T.0) C OUNT { NR ES ) «MI T* KDES-1 
RES(NRES)«FR 
200 CONTINUE 

IF (NIT.LT.MIT.AND.ABS(DG).6T.C0V.AND.ABS(DG) .LT.IO. ) GO TO IbC 
IF ( NOES .GT.O. AND. KDES . EQ. 1. AND. NIT. L T.l) GO TO 150 
IF (NDES.LE.O) GO TO 240 
RDSO*0. 

NDQ»0 

RDQ«0. 

IQ«0 

KQ«0 

DO 210 K=3#MZ 

IF (K.LT.KTEl .0R.K.GT.KTE2) GO TO 210 
CALL VELO (K#K#SV#SM#CP#X#Y#UC#VC#WC) 

11- ITEKK) 

12- ITE2(K) 

CHORD(K ) «X( II )-X(LX ) 

CALL FORCF ( 1 1# 12# X# Y# CP# AL# C HORD ( K ) # XC ( K ) , SC L ( K ) # SCD ( K ) # S CM ( K ) ) 
210 CONTINUE 

CALL TOTFOR ( K TE 1# K TE 2# CHORD# SCL# SCO# SCM# 2 # XC# C L # CDl# CMP# CMR# CM Y# A 
IWING) 

DO 220 I-2#NX 

220 S0{I#2)-3.*<S0(I#3)-S0(I#4))+S0(I#5) 

IF (KDES.LT.NDES) GO TO 140 
GO TO 240 

230 IF (JO.EQ.l) GO TO 10 
JO-1 

GO TO 150 


39 


I 


2^0 RATE«0. 

IF (NRES.GT.l) RATE-{ABS(RES(NRES)/RES(1)) )**( 1 ./( COUNT ( NRES )-CQUN 
1T<1))) 

WRITE (IWRIT>550) 

WRITE {IWRIT#760) RES ( 1 ) / R ES (NR ES ) # COUNT ( NR ES ) > R ATE 
CALL SECOND (T) 

WRITE (IWRIT#770) T 
WRITE (IWRITj670) 

250 CONTINUE 
LX-NX/2+1 
V0L0RG«0. 

DO 260 K-K1#K2 

11- ITEKK) 

X(Il)«XC(K) + .5*SCAL+( A0( 1 1 ) ♦ AO ( 1 1 ) -SO ( I K ) ♦SO ( 1 1# K ) ) 

X(LX)«XC(K) + .5*SCAL*(A0(LX)^A0(LX)-S0(LX,K)tS0(LX>K) ) 
ChORD(K)=X(Il)-X(LX) 

SECDRG(K)»0. 

DO 260 I«2#NX 
260 IDRGPLT( I,K)»-2 
IZDUM-KTE2-1 
S«0. 

DO 270 K»KTE1#IZDUM 
DZ0-.5*(Z(K+1)-Z(K) ) 

270 S-S+DZO* (CH0RD{K+1)+CH0RD(K) ) 

CALL DRAGC (0>0.) 

WRITE (IWRIT#670) 

WRITE (IWRIT,560) VOLDRG 
LX-NX/2+1 

LX0-MIN0(LX+56^ITE2(KT£1)) 

DO 300 I»LX,LXO 
KDUM-0 

DO 280 K«KTE1>KTE2 

280 IF (IDRGPLTdiK) .EU.-2) GO TO 290 
KDUM=KTE2 
GO TO 300 
290 KDUM«K-1 

300 IF (KDUM.GE.KTEl) WRITE (IWRIT#720) ( I DRGP LT ( I^ K ) > K »K TE IMDUM ) 

DO 320 K«3>MZ 

IF (K.LT.KTE1.0R.K.GT.KTE2) 60 TO 320 
I1»ITE1(K) 

12- ITE2(K) 

CALL VELO (K^ K,SV#SM/CP/X,Y,UC^ VC>WC ) 

CHQRD(K)»X( li)-X(LX) 

SECDRG(K) «SECDRG(K) /CHORD (K) 

CALL FORCF (I1»I2#X#Y,CP^AL»CH0RD(K)#XC(K),SCL(K),SCD(K)>SCM(K)) 

IF (KPLOT.GT. l.ANO.K.GT.KTEl) GO TO 310 
WRITE (IWRIT^670) 

WRITE (IWRIT>570) 

WRITE (IWRIT^680) FMACH#YA^AL 
WRITE (IWRITj580) 

310 WRITE (IwRIT^ 680) Z ( K SCL ( K ) >SCD ( K )> SEC DRG ( K )> SC M ( K )> CHORD ( K ) 

IF (KPLOT.LE.l) CALL CPLQT ( I 1> I 2> SM ^ UC # VC > QP RE ( 1^ K ) , A0> SOP RE ( 1> K ) 
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lfSO(l,K)>FMACH) 

320 CONTINUE 

CALL TQTFDR ( KTE 1# KT E 2# CHOR D, SC L>SCD# SCM» XC> C L > C Dl> CMP# CMR> CHY» A 
IWING) 

CD1»CYAW*C01 

CD-CDO+CDl 

VLDl-0. 

IF ( ABS(CDl) .6T.1.E-6) VLDl-CL/CDl 
VLD-0. 

IF (ABS(CD) .GT.l.E-6) VLD=CL/CD 
WRITE (IWRIT#670) 

WRITE (IWRIT#5<J0) 

WRITE (IWRIT#680) FMACH#YA#AL 
WRITE (IWRIT#600) 

WRITE (IWRIT#680) C L# C 0 1# C DO# CD# VLD 1# VL D 
WRITE (IWRIT#610) VOLDRG 
WRITE (IWRIT#620) 

WRITE (IWRIT#680) CMP# CMK# CMY#AWING 
IF (KPLOT.lt. 1) GO TO 330 

CALL THREED ( I P LO T# S V# SM# C P# X # Y # T1 TL E # YA # AL# V L D# C L# CD# CHG R DO# X SC A L 
1#PSCAL#LA6EL#NIT#UC# VC#WC#NF1) 

NFl-11 

IF ( lO.EQ.O) GO TO 230 
330 IF (ISTOP.EO.l) GO TO 380 

IF ( FHALF (NM) .EQ.O. ) GO TO 350 

NX«NX+NX 

NY-NY+NY 

NZ-NZ+NZ 

FDIM=FHALF (2) 

NX1-NX+40-20+FDIM 

NZ1-NZ+2-1+FDIM 

NYl«NY+2-FDIM 

CALL COORD (NX#NY#NZ#KSYM#XTEO#ZTIP#XMAX#ZMAX#SY#SCAL#SCALZ#AX#AY# 
1AZ#A0#A1# A2# A3# B0#Bl#B2#b3#Z#Cl#C2#C3) 

CALL SINGL ( N C#N Z # KS YM# KT E 1# KTE 2# CHOR DC# S W E E P 1# S WE E P2 # S WE £ P# 0 1 HE D1 
1#DIHED2#DIHED#ZS# XLE#YLE#XC#XZ#XZZ# YC# YZ#YZZ# Z»Ci#C2#C3# E1#E2# E3# E 
24#£5# IND#CLO#CLPRE) 

CALL SURF (ND#NE#NC#NX#NZ# IS YM# K SYM# KTE 1# K T E 2# SC AL# YAW# AO# Z#ZS#XC# 
1YC#SL0PT#TRAIL#XS# YS#NP# ITEl# ITE2# IV#SO#ZO# XP# YP#D1#D2»D3# X#Y# IND# 
2XZ#YZ# A1#C1) 

IF (IND.EG.O) GO TO 370 

IF (FDES(l).GT.O. ) CALL SETQS ( NE#N X# QP RE# SO# SOP RE # ITE 1 # I T E2 # K T£1 # 
1KTE2#Z#ZQSTA# A0#PCQl#PCQ2#PCQ3#UC#VC#Q&l#QQ2#QQ3»ClQ4# PCSl#PCS2#DSU 
2RF#NQSTA) 

CALL REFIH (ALFO) 

IREFIN-1 

IF ( 10. EQ.O) GO TO 3A0 

N-Nl 

N1«N2 

N2-N3 

N3-N 

NM-NM+1 
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NIT-0 
GO TO 60 
340 NX-NX/2 
NY-NY/2 
N2-NZ/2 
FDIM=FHALF(1) 

NXl»NX+40-20»FDIM 

N21»NZ+2-l*FDIM 

NY1-NY+2-FDIM 

CALL COORD ( N X» N Y> N2> KS YM# XTEO# 2TI P# XMA X> ZM AX# S Y> SC AL^ SC AL 2# AX# AY» 
lAZ# AO# Al# A2# A3# B0#B1#B2#B3#Z#C1#C2#C3) 

CALL SINGL (N C#NZ# KS YM# KTE 1# KTE 2#CH0RD0# S WE E P 1# S WE E P2# SW E E P# DIME D1 
l#DIriED2#DIHED#ZS#XL£#YLE#XC#XZ#XZZ#YC#YZ#YZZ#Z#Cl#C2#C3#El#£2#E3#E 
24#E5#IN0#CL0#CLPRE) 

CALL SURF (ND#NE#NC#NX#NZ# ISYM#KSYM#KTE 1#KTE2#SCAL> YAW# AO# Z#ZS#XC# 
1YC#SL0PT#TRAIL#XS# YS#NP# ITEl# ITE2# IV»SC#Z0#XP#YP#D1#D2#D3# X#Y# IND# 
2XZ#YZ#A1#C1) 

IF (INO.EQ.O) GO TO 370 
GO TO 230 
350 Kl-KTEl-1 

K2-KTE2+ITE2(KTE2 )-NX/2 

WRITE (8) NX#NY#NZ#NN#K1#K2#NIT 

WRITE (6#630) 

DO 360 K-1#MZ 

WRITE (8) ( (G(I#J#K)#I-1#MX)#J»1#MY) 

360 CONTINUE 

WRITE (8) (E0(K)#K-Ki#K2) 

END FILE 8 
REWIND 3 
GO TO 10 

370 WRITE (IWRIT#67C) 

WRITE (IWRIT#640) 

GO TO 10 
380 CONTINUE 

CALL PLOT (0.#0.#999) 

STOP 0101 
C 

390 FORMAT (28H0NYU INVERSE SWEPT WING CODE) 

400 FORMAT (1X#10A8) 

410 FORMAT (1X#6F10.3) 

420 FORMAT ( 48H0INDICATI0N OF LOCATION OF WING AND VORTEX SHEET#27H IN 

1 COORDINATE PLANE Y - 0 . / 27H0 ( ( I V ( I # K ) # K »K 1# K2 ) » I » 2# N X ) ) 

430 FORMAT ( 49H0C HQRDWI SE CELL DISTRIBUTION IN SQUARE ROOT PLANE#54H A 
IND MAPPED SURFACE COORDINATES AT CENTER LINE AND TIP/15H0 X 

2 #15H ROOT PR0FILE#15H TIP PROFILE ) 

440 FORMAT (15H0 TE LOCATION #15H POWER LAW ) 

450 FORMAT (46H0N0RMAL CELL DISTRIBUTION IN SQUARE ROOT PLANE/15H0 
1 Y ) 

460 FORMAT (15H0 SCALE FACT0R#15H POWER LAW ) 

470 FORMAT ( 4 5H0S P AN W I S E CELL DISTRIBUTION AND SINGULAR LINE/15H0 
1 Z #15H X SING #15H Y SING #15H XZ 

2#15H YZ #15H XZZ #15H YZZ ) 
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480 FORMAT (15H0 TIP LOCATION, 15H POWER LAW ) 

490 FORMAT ( 19H0I TERATIVE S GLUT ION/ 43H0S T R I P WIDTH FOR HORIZONTAL LINE 
1 RELAXATION) 

500 FORMAT (15H0 NX ,15H NY ^15H NZ ) 

510 FORMAT (15H0 MACH NO ,15H YAW ^15H AN6 OF ATTACK) 

520 FORMAT ( lOHOI TE R AT ION, 15H CORRECTION I ,4H J ,4H K ,15H 

1 RESIDUAL I ,4H J , 4H K ,15h SEC LIFT COR ,4H K , lOH 

2 RALF ,10H SONIC PTS) 

530 FORMAT ( lOHOI TERATION, 15H CORRECTION ,4H I ,4H J #4H K ,15H 

1 RESIDUAL ,4H I , 4H J , 4H K , 15H AVERAGE Q ,15H MAXI 

2MUM Q ,4H I ,AH K ,iOH DRAG»10H CL) 

540 FORMAT { lOHOI TE R AT I ON, 1 5H CORRECTION ,4H I , 4H J ,4H K ,15H 

1 RESIDUAL ,4H I ,4H J , 4H K ,10H CL. >10H DRAG ,10 

2H REL FCT 1,10H REL FCT 2,10H BETA ,10H SONIC PTS) 

550 FORMAT (15H0 MAX RESIDAL 1,15H MAX RESIDAL 2,15H WORK ,1 

15H REOUCTN/CYCLE ) 

560 FORMAT ( 13H0WAVL DRAG - ,F9.5,/,23H PRINTOUT OF IDRAG(I,K)) 

570 FORMAT (24HOSECTION CHARACTER1STICS/15H0 MACH NO ,15H Y 

lAW ,15H ANG OF ATTACK) 

580 FORMAT (/13H SPAN S T A T ION, 12 X2HCL, lOX 5HC DO LD, lOX 5HC DN E W, 1 3 X2HC M, 1 0 
1X5HCHQRD) 

590 FORMAT (21H0WING C HAR AC TER I ST IC S / 15ri 0 MACH NO ,15H YAW 

1 ,15H ANG OF ATTACK) 

600 FORMAT (15H0 CL ,15H CD FORM ,15H CD FRICTION ,1 

15H CD ,15H L/D FORM ,15H L/D ) 

610 FORMAT (/2X15H CD WAVE ,/F10.5) 

620 FORMAT (/2X8HCM PITCH, 6X7HCM R0LL,9X6HCM Y A W, 9X SHAW ING ) 

630 FORMAT (1X,14HWRITE ON TAPE8) 

640 FORMAT (24H0BAD DATA, SPLINE FAILURE) 

650 FORMAT (8E10.7) 

660 FORMAT (10A8) 

670 FORMAT (IHl) 

680 FORMAT ( F 12 . 5 , 7F 1 5 . 5 ) 

690 FORMAT (1X,8E15.5) 

700 FORMAT (1H0,10A8) 

710 FORMAT (18,7115) 

720 FORMAT (IX, 3214) 

730 FORMAT ( 1 10, E 15 . 5, 31 4, E15 . 5, 3 14, 5F 10. 5, 1 10 ) 

740 FORMAT ( 1 10, E 1 5 . 5, 3 1 4, E 15 . 5, 3 14, 2E1 5 . 5, 2 1 4, F 10. 4, F 6 . 2 ) 

750 FORMAT ( 1 1 0, E 1 5 . 5 , 3 1 4, E 15 . 5, 3 14 , E15 . 5 , 1 4, F 1 0 . 5, 1 10 ) 

760 FORMAT ( 2E 15. 4, 2F 15 . 4 ) 

770 FORMAT ( 15H0C0MPUTING T IME, FIO. 3,10H SECONDS) 

780 FORMAT ( / 5X3H FNX, IIX 3HFNY , IIX 3HFNZ, 10X5HFPL0T, 9X5HX SC AL, 9X5HPSC AL , 
19X5HFC0NT,10X6HFSWEEP/1X,8£14.5) 

790 FORMAT ( / 4X7HF I T ( NM ) , 8X8HCOVO (NM ) ,8X 7HP10 ( NM ) , 8 X7HP 20 ( NM ) , 8X 7H P30 ( 
INM) , 6X9HBETA0 (NM ),6X9HFHALF (NM) , 6X8HF DE S ( NM ) ) 

800 FORMAT ( /5X5HFMACH, 11X2HY A, 14X2HAL, 12X3HCD0, 12X4HSR EF / IX, 5E15. 5 ) 
END 
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SUBROUTINE GEOM (ND>NC#NP,ZS^XS>YS>XLE^YLE#SLQPT#TRAIL>XP>YP»SWEEP 
11^ SWEEP2>SWEEP>DIHE01>0IHED2>0IHED#XTEC/Ch0RD0>ZTIP»ISYM0>KSYM>CL0 
2) 

C GEOMETRIC DEFINITION OF WING 

DIMENSION XS(ND>1)> YS(ND>1)^ ZS(1)» XLE(D# YLE(D# SLOPT(l)> TRA 
IIL(I)# XP(1)/ YP(1)> NP(1)> CLO(l) 

DIMENSION DESC(IO) 

IREAD-5 

IWRIT«6 

RAD-57,2957795130823 
READ (IREAD^ISO) DESC 

READ ( IRE AD> l<fO) ZSYM# FNC > SwE E P 1> S WE E P2 # S WE E P # D I HE D D IHE D2> D IHE 0 

WRITE ( IWRIT> 190) ZS YM> FNC# SWEE Pl> S WE E P 2# S W E E P> D IHED1> D IHE D2> DIHE D 

IF (FNC.LT.3. ) RETURN 

KSYM-ZSYM 

NC-FNC 

SWEEPl-SWEEPl/RAD 

SWEEP2-SWEEP2/RAD 

SWEEP-SWEEP/RAD 

DIHEDl-DIHEDl/RAD 

DIHED2-DIHED2/RAD 

DIHED-DIHED/RAD 

ISYMO-1 

XTEO-0, 

ChORDO-0. 

K-1 

10 READ {IREAD/150) DESC 

READ (IREADil^O) ZS ( K ) > XL # YL# CHORD# TH ICK # A L> FSEC > C LO ( K ) 

WRITE (IWRIT#200) Z S ( K ) # X L # Y L # CHORD# TH I CK# A L# F S E C 
ALPHA-AL/RAD 

IF (K.GT.l. AND.FSEC .EQ.O. ) GO TO 80 
READ (IREAD#150) DESC 

READ (IREAD#140) YS YM# F NU# FN L # SNGOPT# ZE ND# S NGRA T 

WRITE (IWRIT#210) YSYM#FNU#FNL 

NU-FNU 

NL-FNL 

N-NU+NL-1 

READ (IREAD#150) DESC 
READ (IREAD#140) TRL# SLT# XSING# YSING 
WRITE (IWRIT#220) TR L# S L T# XS ING# YS I NG 
READ (IREAD#150) DESC 
WRITE (IWRIT#230) 

DO 20 I-NL#N 

READ (IREAD#1A0) XP(I)#YP(I) 

20 WRITE (IWRIT#180) XP(I)#YP(I) 

L-NL+1 

IF { YSYM .GT.O . ) GO TO 40 
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READ (IREAD#150) DESC 
WRITE (IWRIT>240) 

DO 3.0 I«1»NL 

READ (IREAO/140) VAL,DUM 
WRITE (IWRIT#180) VAL^DUM 
J«L-I 
XP( J ) *VAL 
30 YP(J)»DUM 
GO TO 60 
^0 J-L 

DO 50 I»NL#N 
J-J-1 

XP(J)»XP{I) 

50 YP(J)»- YP(I) 

60 WRITE (IWRIT/160) 

WRITE (IWRIT^IIO) ZS(K) 

WRITE (IWRIT#170) TRL#SLT# XSING> YSING 
WRITE (IWRIT,120) 

DO 70 I*1>N 

70 WRITE (IWRIT>170) XP(1),YP(I) 

80 CONTINUE 

SCALE»CH0RD/(XP(1)-XP(NL)) 

XLE(K) *XL+( XSING-XP (NL ) )*THICK*SCALE 
YLE (K) -YL+ ( YS ING-YP (NL ) ) ♦ T HIC K* SC AL E 
XX»XP(NL)+(XS IN G-XP(NL))* THICK 
YY»YP(NL)+ (YS ING-YP (NL ) )*THICK 
CA»COS( ALPHA) 

SA«SIN( ALPHA) 

DO 90 I»1#N 

XS(I#K) -SCALE *( (XP( I)-XX)*CA+THICK»(YP( I)-YY)*SA) 

90 YS( I>K )» SC ALE* (THICK* (YP( I)-YY)*CA-(XP(I)-XX)*SA) 

SLOPT(K ) -THICK*SLT-T AN (ALPHA) 

TRAIL (K)»THICK*TRL/R AD 
NP(K)«N 

XTE0-AMAX1(XTE0/XS(1^K)) 

CHORDO-AMAXK C HO RDO^ CHORD ) 

IF ( YSYM.LE.O.. OR. ALPHA. NE.O. ) ISYMO-0 
WRITE (IWRIT,130) ZS(K) 

WRITE (IWRIT*170) XL# YL ^ CHORD/ T H ICK> AL 
K-K + 1 

IF (K.LE.NC) GO TO 10 
Z0-.5*( ZS( 1) + ZS(NC ) ) 

IF (KSYM.NE.O) ZO-ZS(l) 

DO 100 K-l/NC 
100 ZS(K)-ZS(K)-ZO 
ZTIP-2S(NC ) 

RETURN 

C 

110 FORMAT (16H0PR0FILE AT Z - /F10.5/15H0 TE ANGLE /15H TE SL 

lOPE /15H X SING /15H Y SING ) 

120 FORMAT (15H0 X /15H Y ) 

130 FORMAT (27H0SECTI0N DEFINITION AT Z - /F10.5/15H0 XLh »15 
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IH 

2PHA 

140 FORMAT 
150 FORMAT 
160 FORMAT 
170 FORMAT 


YLfc 


15H 


CHORD 


» 15HTHICKNESS RATIDilSH 


AL 


) 


(8F10.6) 

(10A8) 

( IHl) 

( F12.4,7F15.4) 

180 FORMAT (8E15.5) 

190 FORMAT ( /5X4HZSYM> 12X3HFNC # 10X6HSWE EP 1> 9X6HS WEEP 2> 9X5HS WEE Pj 10X6HD 
IIHEDI/ 9X6HDIHED2> lOX 5H D IHE 0/ 1 Xj 8E15. 5 ) 

200 FORMAT ( / 5X5HZS ( K ) ^ 1 2X2HX L# 13 X2HYL# 11 X5HCH0RD/ lOX 5HTH ICK> 1 2X2HAL> 1 
12X4HFSEC/1X,7E15. 5) 

210 FORMAT ( / 6X4H YS YM# 11 X3HFNU> 12X3HFNL / IX/ 3E 1 5 . 5 ) 

(/6X3HTRL>12X3HSLT/llX5HXSING/lCX5hYSING/lX/4E15.5) 
(/5X5HXP(I)/10X5HYP(D) 

(/6X3HVAL/12X3HDUM) 


220 FORMAT 
230 FORMAT 
240 FORMAT 
END 


SUBROUTINE COORD (NX/NY/NZ/KSYM/XTEO/ZTIP/XMAX/ ZMAX/ S Y/ SC AL/ SCALZ# 
lAX/ AY/ AZ/ AO/ A1/A2/A3/ BO/ 81/ B2/B3/Z/C1/C2/C3) 

C SETS UP STRETCHED PARABOLIC AND SPANWISE COORDINATES 

DIMENSION A0(1)/ Al(l)/ A2(l)/ A3(l)/ B0(1)/ Bid)/ 82(1)/ B3(l)/ 
1Z(1)/ Cl(l)/ C2(l)/ C3(l) 

COMMON /DIM/ NXl/NYl/NZl/ FDIM 

DX-2./NX 

DY-l./NY 

KY-NY+1 

DZ-2. /NZ 

DX«2. /NXl 

DY-l./NYl 

DZ-2./NZ1 

KY1»NY1+1 

ZO-l.-DZ 

Kl«2 

K2»NZ 

IF (KSYM.EU.O) GO TO 10 

DZ-l./NZ 

DZ-l./NZl 

Z0»0. 

Kl«3 
K2-NZ+2 
10 AX«,5 
AY-. 5 
AZ». 5 
BX«0, 

BZ«0. 
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XMAX-.625 
ZMAX»8. /lA. 

SY«. 5 

SCAL«XTEO/ ( ,50001+XMAX+XMAX) 

SCALZ»ZTIP/ (l.OOOOOl^ZMAX) 

V2»(DX/DY)**2 
Wl-SCAL/SCALZ 
W2»(W1*DX/DZ)**2 
S73-SQRT(73. ) 

BBX— BX+SaRT(3.*(7.+S73) ) /( (l.+S73)*XMAX**3) 

ABX»1.-BBX*SQRT( ( 7 . +S 73 ) / 12 . ) *XM AX*+ 3 
CBX« ( 19.+S73) ♦XMAX+XMAX/IZ. 

ABBX-ABX+BBX* (3 .♦CBX-^.*XMAX*XHAX)+XMAX+XMAX/SQRT(CBX-XMAX*XMAX) 

MX-NX+1 

DO 40 I-l/MX 

DD»( I-l)*DX-l.+(NXl-NX)+DX/2. 

B-1. 

IF {ABS(OD) .GT.XMAX) GU TO 20 

A-CBX-DD+DD 

AS-SORT(A) 

C»ABX+AS+BBX* ( 3 . * C BX-4 . * DD+ DD ) * DD*OD 

DU»ABX*DD+8BX+AS+DD**3 

Dl-AS/C 

D2-SBX*(CBX*(-6.+CBX+19.4DD+DD)“12.*DD«+4)*D0/(A*C) 

GO TO 30 

20 IF (DO.LT.O.) B»-l. 

A«l.-((DD-B*XMAX)/(l.-XMAX))+^2 

C-A++AX 

D»{AX+AX-1. )♦ ( l.-A ) 

D0*8*XMAX+ABBX*( OD-B+XMAX ) /C 
01»A*C/( (l.+D)+ABBX) 

D2»-{ AX+AX)*{ DD-B+XMAX)*( 3.+D )/ ( ( 1 . +D ) ♦ A« ( 1 . -XM A X ) **2 ) 

30 AO(I)-DO 

Al( I)».3*D1/DX 
A2(l )»D1*D1 
A3( I ) ». 5+DXtD2 
40 CONTINUE 
JS«4-FDIM 
DO 50 JJ*JS^KY1 
J»JJ-( 2-FDIM) 

DD»(KY1-JJ )+DY 

A=1.-DD*D0 

C«A**AY 

D«{ AY+AY-1, )♦ (l.-A) 

Dl-A+C/{ (l.+D)*SY) 

B0( J ) «SY*DD/C 
BKJ )«.5*D1/DY 
B2( J ) *D1*D1*V2 

50 B3(J )— AY+0D*DY+(3.+D)/( (l.+D)*A) 

BBZ»-BZ*SQRT( 3.* ( 7.+S73 ) ) / ( ( 1 . + S73 ) ♦ ZM A X*+ 3 ) 

ABZ=1.-BBZ+SQRT( ( 7. +S73 ) / 12 . ) ♦ZMAX** 3 
C BZ» ( 19. +S73)*ZMAX*ZM AX/12. 
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ABBZ«ABZ + BbZ* (3.»CBZ-<f.+ZMAX*ZMAX) + ZMAX + ZMAX/SQRT(CBZ-ZMAX*ZMAX ) 

DO 80 K»2^K2 

DD*(K-K1)«DZ-Z0 

B*1 • 

IF ( ABS ( DO ) .GT.ZMAX ) GO TO 60 

A-CBZ-DO+DD 

AS*SQRT( A) 

C»ABZ*AS+BBZ* ( 3 . +C BZ-^ . *DD+ DO )*OD*DD 

D0-ABZ*DD+BBZ+AS*D0**3 

Dl-AS/C 

D2-BBZ+ (CBZ*(-6.*CBZ + 19.*DD*D0)-12. + DD* + A)*DD/ ( A*C ) 

GO TO 70 

60 IF (DO.LT.O.) B«-l. 

A-l.-{ ( DO-B + ZMAX) / (l.-ZMAX ) )**2 
C«=A**AZ 

D»(AZ + AZ-1. )>M1.-A} 

D0*B>('ZMAX + ABBZ*(DD-B+ZMAX )/C 
Dl«A*C/( (l.+D)«ABBZ) 

D2«-(AZ+AZ)*(DD-B^ZMAX)*(3.+D)/((1.+D)*A*(1.-ZMAX)*^2) 

70 Z(K) =SCALZ’>DO 

C1(K)«.5*D1+W1/DZ 
C2(K)-D1*D1*W2 
C3(K)«.5*DZ*D2 
80 CONTINUE 
RETURN 
END 


SUBROUTINE SINGL { NC > N Z^ KS YM# KTE 1# K TE 2# C HOR DO/ S WEEP S WE E P 2, S kJE EP , 
1DIHED1»DIHED2#DIHED/ZS>XLE>YLE#XC*XZ>XZZ#YC^YZ#YZZ>Z>C1#C2>C3#E1#E 
22#E3^E4^ E5# IND^CLOi CLPRE) 

C GENERATES SINGULAR LINE FOR SQUARE ROOT TRANSFORMATION 

DIMENSION ZS(1)^ XLE(l)f YLE(D# XC ( 1 ) # X2(D# XZZ(1)» YC(i)> YZd 
1)> YZZ(1)> Z(l)> Ci(l)j C2(l)> C3(D# El(l), £2(1)> E3(D# EA(D* 
2E5(l)^ CLO(D# CLPRE(l) 

DO 10 K»1,NC 
EA(K)»0. 

10 E5(K)»0. 

Kl-2 

K2-NZ 

IF (KSYM.EQ.O) GO TO 20 

Kl«3 

K2-NZ+2 

KTEl-3 

20 DO 3C K=K1^K2 

IF (Z(K).LT.ZS(D) KTEl-K + 1 
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IF ( Z(K) .LE.ZStNC ) ) KTE2-K 
30 CONTINUE 
B-CHQRDO 
S1»TAN(SWEEP1 ) 

S2»TAN(SWEEP2) 

T1«TAN(0IH£D1 ) 

T2«TAN(DIHE02) 

CALL SPLIF (l^NC^ ZSfXLE^ E1#E2/E3^1^ S1»1#S2^0>0.^ IND) 
CALL INTPL (KTE1#KTE2> Z^XC#1>NC#ZS>XLE> E1^E2>E3»0) 
CALL INTPL lKTEl,KTEZ,Zf)(Z,l,HC,ZS,BlfBZ,E3,t^,0) 

CALL INTPL (KTE1>KTE2>£»XZZ>1>NC>ZS» E2>E3>£<»>E5>C) 
CALL SPLIF <1>NC#ZS#YLE#E1,E2#E3>1^T1#1>T2>0>0.> IND) 
CALL INTPL (KTE1^KTE2^Z^YC>1,NC^ZS,YLE#E1#E2>E3,0) 
CALL INTPL (KTEl>KTE2iZ#YZ#l#NC/ZS>ei#E2>E3>E^#0) 

CALL INTPL (KTEl,KTE2*Z»YZZ>l,NC>ZS»E2,E3,EA,Eb,0) 
CALL SPLIF C 1 ^ NC > ZS ^ C LO » E 1 # E 2> E 3# 3> 0 • j 3> 0 • # 0> 0 • / I NO ) 
CALL INTPL (KTE1>KT£2^Z/CLPRE^1^NC,ZS>CL0#E1>E2^ E3>0) 
S«B*TAN ( SWEEP ) 

Sl-B^Sl 

SZ-BX'SZ 

T-B*TAN(OIHED) 

Tl-B+Tl 

TZ-B^TZ 

XC(2) «3.+( XC( 3)-XC (4) )+XC (b) 
YC(Z)-3.+(YC(3)“YC(4))+YC(b) 

IF (KSYM.NE.O) GO TO 50 
N-KTEl-1 
DO 40 K-K1,N 
ZZ-(Z(K)-Z(KTE1) )/B 
A-EXP(ZZ) 

XC(K) ■XC(KTE1 )+S+ZZ-( Sl-S )+(l.-A) 

YC(K )-YC (KTEl )+T*ZZ-(Tl-T)*<l.-A) 

XZ(K)»(S+(S1-S)*A)/B 

YZ(K)-(T+(T1“T)+A)/B 

XZZ(K)»(S1-S)*A/(B*B) 

40 YZZ(K) «( Tl-T) *A/ (B + B ) 

50 N-KTE2+1 

DO 60 K«N^K2 

ZZ-( Z(K)-Z(KTE2) ) /B 

A-EXP(-ZZ) 

XC(K)»XC(KT£2 )+S*ZZ+(S2-S )*(!.- A) 

YC (K)“YC (KTE2 )+T*ZZ+(T2-T)*(l.-A) 

XZ(K)«(S+(SZ-S»*A)/B 
YZ(K)«(T+(TZ-T)*A)/B 
XZZ(K) *-(S2-S )*A/ (B + B ) 

60 YZZ{K)«-(T2-T)*A/(B*B) 

RETURN 

END 
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SUBROUTINE SURF { NO# NE# NC# NX# N2 > ISYM# KS YM^ K7 E 1# KTE 2, SC AL, Y AW^ AO* Z, 
1ZS#XC# YC#SLQPT#TRAIL#XS#YS#NP#ITE1#ITE2#IV#SO#ZO#XP#YP#D1#02#03#X# 
2Y# IND#XZ# YZ#A1#C1) 

C INTERPOLATES MAPPED WING SURFACE AT MESH POINTS 
C INTERPOLATION IS LINEAR IN PHYSICAL PLANE 

DIMENSION SO{NE#D# XS(ND#D# YS(ND#D# ZS(D# SLOPT(l)# TRAIL(l)# 
1 XC(D# YC(D# AO(D# Z(D# ZO(D# X(D# Y(D# XP(D# YP(D# 01(1) 
2# D2(D# 03(1)# 1V(NE#D# NP(1># ITEKD# IT£2(i)# XZ(D# YZ(D# A 
31(1)# CKl) 

COMMON /DIM/ NX1#NY1#NZ1#F0IM 
PI-3.14159265268979 
TYAW-TAN( YAW) 

S1-. 5+SCAL 

DX-2./NX1 

LX-NX/2+i 

MX-NX+1 

MZ-NZ+3 

IVO-l-ISYM-ISYM-ISYM 
IV1--1-ISYM 
DO 1C K»1#MZ 
ITEKK) =MX 
ITE2(K) -MX 
DO 10 I-1#MX 
IV ( I#K) =-2 
10 S0(I#K)«0. 

K-KTEl 

K2-1 

20 K2-K2+1 
K1-K2-1 
R2-1. 

IF ( ZS(K2)-Z(K) ) 20#40#30 
30 R2-{Z(K)-ZS(K1) )/(ZS(K2)-ZS(Kl) ) 

40 R1-1.-R2 

C-Rl*XS(i#Kl)+R2*XS(l#K2) 

CC-SQRT( (C+C) /SCAD 
DO 50 I»2#NX 

IF ( ( A0( I ) + .5*DX ) .LT.-CC ) 11 = 1+1 
IF ((A0(I)-.5*DX).LT.CC) 12-1 
50 CONTINUE 
ITEKK) -II 
ITE2(K) =12 
CC»AO( 12) /CC 

Z0(K)-Z(K)-TYAW*(XC(K) + S1*A0( I2)*A0( 12) ) 

KK»K1 

P-Rl 

60 N-NP(KK) 

Q-SQRT(XS(1#KK)/C)/CC 
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DO 70 I»2>NX 
70 X{I)»Q*AO(I) 

ANGL-PI+PI 

U-1. 

V-0. 

DO 90 I»1>N 

R*SQRT(XS(I>KK)**2+YS(I/KK)*+2) 

IF (R.EQ.O.) GO TO 80 

AN6L»ANGL+ATAN2( (U*YS( UKK)-V*XS( I>KK) ) > ( XS ( I > KK ) + V* YS ( I# KK ) ) ) 
U«XS( I^KK) 

V«YS(I,KK) 

R-SQRT( (R+R)/SCAL) 

XP( I)-R*C0S{.5*ANGL) 

YPd ) »R*SIN( . 5+ANGL ) 

GO TO 90 
80 ANGL«PI 
U»-l . 

V»0. 

XP(I)»0. 

YP(I)-0. 

90 CONTINUE 

ANGL = ATAN( SLOPKKK ) ) 

ANGL1-ATAN(YS (1>KK)/XS(1/KK) ) 

ANGL2-ATAN(YS(N#KK)/XS(N>KK)) 

ANGL 1-ANGL-. 5* (ANGL 1-TRAIL (KK) ) 

ANGL2*ANGL-.5*(ANGL2+TRAIL(KK)) 

T1«TAN( ANGLl) 

T2«TAN( ANGL2) 

CALL SPLIF (l#N>XPj YPj01iD2jD3^1>Tl#l>T2*0>0.#IND) 

CALL INTPL (I1#I2#X^Y#1>N»XP#YP^D1#D2>D3»0) 

X1-.25+XS ( 1/KK ) 

A«SLOPT(KK)+( XS( 1>KK)-X1) 

B-1./(XS(1#KK )-Xl) 

ANGL-PI+PI 

U»l. 

V«0. 

M»I1-1 

DO 100 I*2>M 
XX«. 5*SCAL+X( I)»*2 
D«B* ( XX-Xl ) 

YY»YS(1/KK)+A*AL0G(D)/D 
R»SQRT( XX++2+YY**2) 

ANGL»ANGL + ATAN2( ( U+YY-V + XX ) > CU*XX+V*YY) ) 

U«XX 

V-YY 

R»SQRT( ( R + R)/SCAL ) 

100 Y(I)*R*SIN(.5+ANGL) 

A-SLOPT(KK )♦( XS(N>KK}-X1) 

B«1./(XS(N#KK )-Xl) 

ANGL-0. 

U»l. 

V»0. 
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II II II Mil 


M»I2+1 

DO 110 I»M^NX 
XX«.5+SCAL*X( l)*+2 
D»B* ( XX-Xl ) 

YY*rS(N# KK)+A=t‘ALOGlD)/D 
R-SQRT( XX%*2+YY**2) 

ANGL- ANGL+ATAN2( ( U* YY-V+XX ) # ( U+XX+V+ YY ) ) 

U«XX 

V«YY 

R-SQRK (R+R)/SCAL) 

110 Yd )=R*SIN( .5*ANGL ) 

Q-P*Q*CC*CC 
DO 120 I«2>NX 

120 S0(I#K)-SC(I#K)+Q*Y(1 ) 

IF (KK.EQ.K2) GO TO 130 

KK-K2 

P«R2 

GO TO 60 

130 DO lAO I«=I1>I2 
140 IV(I>K)«2 
M-Il-1 

00 150 I*2>M 

ZZ«Z(K)-TYAW*(XC(K)+S1*A0(I)+A0(I)) 

IF (ZZ.GE.ZO(KTEl) ) IV(I#K)«IV0 
150 CONTINUE 
M»I2+1 

DO 160 I«M#NX 

Z2-Z(K)-TYAW* (XC(K)+S1*A0( I)*AO(I) ) 

IF (ZZ.GE.ZO(KTEl) ) IV(I>K)«IVO 
160 CONTINUE 
K2“K2-1 
K-K + 1 

IF <K.LE.KTE2) GO TO 20 

Kl = 2 

K2-NZ 

IF (KSYM.EQ.O) GO TO 170 

Kl»3 

K2»NZ+2 

170 DO 160 I«2#NX 

ZZ»Z(K)-TYAW* (XC (K)+S14A0( I)*AO( I) ) 

IF (ZZ.LE.ZS(NC) .AND.ZZ.GE.ZO(KTEl) ) IV(I#K)»IVO 
180 CONTINUE 
K-K + 1 

IF (K.LE.K2) GO TO 170 
N-KTE2 

IF ( YAW.LE.O. ) GO TO 200 
I0-ITEKKTE2) +1 
DO 190 I-IO»LX 
N-N + 1 

190 Z0(N>-Z(KT£2)-TYAW*(XC(KTE2)+S1*A0< I)*AO(I) ) 

200 I-ITEKKTEl) 

Z0(KT£1-1)«Z(KTE1-1)-TYAW*(XC (KT E 1-1 ) +S !♦ AO ( I )♦ AO ( I ) ) 
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Z0(N+1)»Z{KTE2+1) 

DO 220 K«K1>K2 
00 210 I-2/NX 

IF (IV(l#K).GT.O) GO TO 210 

IF (IV(I + l>K-H).GT.O.OR.IV(I-l»K + l).GT.O) IV(I#K)«IV1 
IF ( IV( I+1*K-1) .6T.0.QR.IV(I-1>K-1) .GT.O) IV(I>K)-IV1 
210 CONTINUE 

220 IF (S0(LX>K).LT.l.E-05) IV(LX>K)>0 
IF (KSYM.EQ.O) RETURN 
00 230 I«2>NX 

230 S0(I>2)»3.*(S0(I#3)-S0( I»4))+S0(I,5) 

RETURN 

END 


SUBROUTINE ESTIM (ALFO) 

C INITIAL ESTIMATE OF REDUCED POTENTIAL 

COMMON G(129^ 12> 15 ) > S0( 129# 15 )> E0( 131 ) / ZO ( 131 )> I V ( 129# 15 ) ^ I TEl ( 15 ) 
1#ITE2(15)#A0( 129)# A1(129)#A2(129)#A3(129)#B0(12)#B1<12)#B2(12)#B3( 
212)#Z(15)#C1(15)#C2(15)#C3(15)#XC(15)#XZ(15)#XZZ(15)#YC(15)#YZ(15) 
3# YZZ(15)#NX#NY#NZ#KTE1#KTE2# IS YM# KS YM, S C AL # SC AL Z# Y A W# C Y A W, S Y A W # ALP 
4HA#CA#SA#FMACH#N1#N2#N3#I0#ND£S#TSTEP#EPS1#QPRE(129#15)#S0PRE(129# 
515)#NQSTA#ZQSTA(15)#PCQ1(15)#PCQ2(15)#PCQ3(15)#QQ1(15)#QQ2(15)#QQ3 
6(15)#QQ^(15)> PCS1(15)#PCS2(15)#DSURF(15)>RDQ#R0S0#F00#F01#F10# Fll* 
7N0Q# IQ#KQ# AWING# VOL DRG# lORGP LT( 129# 15 ) # SECDRG ( 1 5 ) 

DIMENSION ALFO(l) 

KY-NY+1 
MZ-NZ+3 
DO 10 I»l#129 
DO 10 J«l#12 
DO 10 K»l#15 
10 G(I#J#K)>0. 

K»1 

DO 30 K-1#MZ 
DO 20 1-2#NX 
G(I#KY+1#K)«0. 

IF ( IV< I#K) .LT.2) 60 TO 20 
DSI-SO( I+1#K)-S0( I-1#K) 

DSK-S0( I#K+1)-S0( I#K“1) 

SX-A1(I)*0SI 

SZ»C1(K)*DSK 

FH-A0( I)*A0(I ) + S0(I#K) + S0( I#K) 

H«1./FH 

AZ— A0{I)*XZ(K)-S0( I#K)*YZ(K) 

BZ*-A0( I)*Y2(K)+S0( I#K)*XZ(K) 

HZ«AZ*SX-BZ+FH*SZ 
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FYY«1.+SX*SX+H*HZ*HZ 
FXY-SX+H+AZ+HZ 
VSA+AO( I)-CA+SO( I^K) 

U»CA*AO(I)+SA*SO(I>K) 

W-SYAW+CA*XZ(K)+SA*YZ(K) 

G( I>KY+1»K)«G (I#KY-l#K) + ( V*a.-H*8Z*HZ)-U^FXY-W*HZ) / (FYY + BKKY) ) 
20 CONTINUE 
30 CONTINUE 
Kl-KTEl-1 

K2«KTE2+ITE2(KTE2)“NX/2 
DO 40 K»K1>K2 
ALFO(K) >0. 

40 E0(K)»0, 

I0»1 

RETURN 

END 


SUBROUTINE MIXFLQ 

C SOLUTION OF EQUATIONS FOR MIXED SUBSONIC AND SUPERSONIC FLOW 
C USING ROTATED DIFFERENCE SCHEME 

COMMON G(129, 12^ 15 ) # S0( 129, 15 ) ^ E 0( 131 ) > ZO ( 131 ) » I V ( 129, 15 ) , ITEl ( 15 ) 
1,ITE2(15), A0( 129),A1(129),A2(129),A3(129), B0(12),B1(12),B2(12),B3( 
212),Z(15),C1( 15),C2(15),C3(15),XC(15),XZ(15),XZZ(15),YC(15),YZ(15) 
3, YZZ( 15),NX,NY#NZ,KTE1,KTE2, ISYM,KSYM,SCAL, SC AL Z , Y AW, C YA W, S YAW, ALP 
4HA,CA,SA,FMACH,N1,N2,N3,IQ,NDES,TSTEP,EPS1,QPRE(129,15),S0PRE(129, 
515),NQSTA,ZQSTA(15),PCQ1(15),PCQ2(15),PCQ3( 15),QQ1( 15),QQ2(15),QQ3 
6(15),QQ4(15),PCS1(15),PCS2(15),DSURF( 15),R0Q,RDS0,F00,F01,F10,F11, 
7NDQ, IQ,KQ, AW I NG, VOL DRG, IDRGPLT( 129, 15 ) , SEC DRG ( 1 5 ) 

COMMON /FLO/ STR I P, PI, P 2, P3, BET A, FR, I R, J R, K R, DG, IG, JG, KG, NS, PS WEE P 
COMMON /SWP/ DXYZ(129),GK1(129,15),GK2(129,15),SX(129),SZ(129),SXX 
1(129),SXZ(129),SZZ(129),R0(129),R1(129),C(129),D(129),G10(15),G20( 
215),G30(15),G40{15),G1{15),G2(15),I1, I2,K, L ,N0, L X, M X,K Y, MY, T 1, AAO, 
3Q1,Q2,TYAW,S1 

COMMON /DIM/ NX1,NY1,NZ1,FDIM 

BETX-.Ol 

BETY-.15 

BETZ-.l 

BSCAL«1./(1.+FDIM) 

BSCAL1«1./(2.*<1.+FDIM) ) 

LX-NX/2+1 

MX-NX+1 

KY-NY+1 
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HY»NY+2 

TYAWSYAW/CYAW 

S1».5*SCAL 

DX-2./NX1 

T1»DX*0X 

AA0-l./FMACH**2+.2 
Ql-2. /PI 
Q2»l. /P2 
FR»0. 

IR-0 

JR»0 

KR»0 

DG-0. 

IG«0 
J G»0 
KG-0 
NS«0 
Kl-2 

IF ( FMACH.GE. 1. ) Kl-3 
K2»NZ 

IF (KSYM.EQ.O) GO TO iO 

Kl-3 

K2-NZ+2 

10 F«A6S( . 5+STRIP+NX) 

L*F 

IF (L.EQ.NX/2) L-L-1 

I1«LX-L 

I2-LX+L 

IF (L.EQ.O) I2»LX-1 
DO 20 J»1,MY 
DO 20 I-l^MX 
GKK I>J )»G( I» J#l) 

20 GK2(I>J)*G( li Jjl) 

K«2 

L-2 

N0»KTE1-1 

IF (K.EQ.Kl) GO TO 90 
IF (KSYM.EQ.O) GO TO 80 
I»LX 

DSI-S0( I+l#3)-S0( I-l#3) 
DSK-S0(I#<*)“S0(I#2) 

SX(I )-Al( I )*DSI 
SZ( I)-C1(3)*DSK 
R»1 . 0 

DO 30 J«2>KY 
YP»BO(J )+S0(I#3) 

IF (J.EQ.KY) R«AMIN0( i>IV ( I>K ) ) 
H-R/(1.-R+YP*YP) 

AZ»-YP*YZ ( 3 ) 

BZ«YP*XZ(3) 

A-H+AZ+AK I) 

B«(H+(BZ-AZ*SX(I))-SZ(I))*B1(J) 



DGI»6( I+l# J>3 )-G( I-l# J>3) 

DGJ*G(I j J+l#3)-G( I> J-l#3) 

G( I> J^2)«G(I> A)+(A+DGI~B+DGJ)/C1C3) 

GKK 1,J )»G( I, J#2) 

G(I>J#l)«3.+(G(I#J#2)-Ga,J,3))+G(I,J,A) 

GK2( I# J)«G(I, J^l) 

30 CDNTINUfc 
J»KY+1 

G(I>J>2)=G(I^ A) + (A*DGI-B*DGJ) /Cl( 3) 

GKK I#J)«G(I, J>2) 

G(I,J#l)«3.+(G(l#J>2)-G(I,J/3))+G(I>J>A) 

GK2(I/J)=G(I/ J^l) 

M-NX/2-1 
DO 7C 

GO TO 50 
40 I-LX+II 

50 DSI»SO( I+l#3)-S0( 1-1^3) 

OSK=SO( I>4 )-S0(l#2) 

SX( I ) «A1 ( I )*DSI 
SZ(I)«C1(3)*DSK 
DO 60 J*2#KY 
YP»B0(J)+S0(I^3) 

H« 1. /( AO ( I )*A0( I )+YP*YP) 

AZ»-AO(I )*XZ( 3)-YP+YZ(3) 

BZ — A0( I)*YZ(3)+YP*X2(3) 

S-SIGN ( 1 AZ) 

A-H*ABS( AZ)*A1(I ) 

B«(H=«'( BZ-AZ*SX( I ) )-SZ{I))*Bl(J) 

IP»I+IFIX(S) 

IM*I-IFIX(S) 

OGI*G( I,J^4)-G( IM> J,4) 

0GJ«G(IfJ+l»3)-G(I^J-l>3) 

G(I>J»2) = (Cl(3)*G(I>J#‘t) + A*(G{IP,J,2)+DGI)-B+DGJ)/{Cl(3)+A) 
GKK K J ) «G( K J,2 ) 

G(I#J>1)»3.*(G(I>J>2)-G(KJ>3))+G(I#J>4) 

60 GK2(KJ)=G(I> J>1) 

J-KY+1 

G(KJi2)»(Cl(3)*G(KJ#4) + A*(G(IP#J,2)+DGI)-B+DGJ)/(Cl(3)+A) 
GK1(KJ)=G(KJ,2) 

IF ( I.LT.LX) GO TO 40 
70 CONTINUE 
80 KK-K+1 
K3-K2+1 

90 DU 150 K*KK#K2 
DO ICO J«1>MY 
G10( J )-G( I2#J #K ) 

G20( J)»G(I2-1#J/K) 

G30(J)-G(I1>J>K) 

100 G40( J )«G( I1+1#J>K) 

DO 110 I»2*NX 

DSI»S0( I+1>K)-S0( I-1,K) 
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DSK«SO( I,K+1)-S0( I#K-1) 

DSII«S0{I+1#K)-S0(I»K)-S0(I»K)+S0(I-1>K)+A3(I)*DSI 
DSKK»SO( I/K+1)-S0(I#K)-S0{I#K)+S0(I#K-1)+C3(K )+DSK 
DSIK*SO( I+1>K+1)-S0( I-1>K+1)-SG(I+1>K-1)+S0{ I-l/K-1 ) 
SX(I)-Al(i)+DSl 
SZ(I )=C1(K)*DSK 
SXX( I)«A2(I)^DSII 
SZZ( I)-C2(K)+DSKK 
110 SXZ(U=T1*A1( I)*CUK>*DSIK 
L«K 

IP (I2.LE.il) GO TO 130 
IF (FSWEEP.LT.O. ) GO TO 120 
CALL YSWEEP 
GO TO 130 
120 CALL VYSWEEP 
130 CONTINUE 

IF (K.NE.KTE2.0R.YAw.Le.0. ) GO TO 150 

I0-ITE1(K)+1 

DO lAO I=10#LX 

M-NX+2-I 

E«G(M>KY#K)-G(I#KY>K) 

NC-NO+1 

lAO E0(N0)»E0(N0)+P3*(E-E0(N0) ) 

150 CONTINUE 

C BOUNDARY CONDITION AT INFINITY REPLACED BY MIXED DIRlChLET 
C AND NEUMANN CONDITION AT CONTROL SURFACE 
DO 160 I-1#MX 
DO 160 J-1,MY 

160 G(Ij J^K 3)* ( 1. -BETZ/BSC AL1)*G(I>J#K2) 

DO 170 J«1>MY 

G{I2+1#J#2)-(1.-BETX/BSCAL)*G(I2#J>2) 

170 G(I1-1>J#2)=(1.-6ETX/BSCAL)*G(I1^J#2) 

DO 180 I»1^MX 

180 G(I,J1-1#2)»(1.-BETY/BSCAL1)*G(I,J1,2) 

FR«1.2*FR/ AAG 

RETURN 

END 


SUBROUTINE YSWEEP 

C THE FINITE DIFFERENCE EQUATIONS FOR G ARE SOLVED BY ROW RELAXATION 

C MOST OF THE COMPUTING TIME IS SPENT IN THIS ROUTINE 

COMMON G(129, 12#15)jS0(129,15)»E0(131)»Z0(l31)>IV(12‘y#15)>IT£l(15) 
1>ITE2(15)> A0(129),A1(129)#A2(129),A3(129),B0(12)^B1(12),B2(12)>B3( 
212)»Z{15),C1(15)>C2(15)>C3(15)^XC(15)>XZ(15)>XZZ(15)>YC(15)>YZ(15) 
3#YZZ(15)/NX#NY>NZ»KTE1#KTE2#ISYM,KSYM>SCAL#SCALZ>YAW»CYAW>SYAW,ALP 
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4HA>CA>SA» FfiACH#Nl,N2#N3#IQ>NDES» TSTEP#EPS1#QPRE(129,15)* S0PR£(129> 
515)/NQSTA#ZQSTA(15)#PCQ1(15)#PCQ2(15 )>PCQ3( ) » QQl ( ) » QQ2 ( 15 ) > QQ3 
6(15),QQA(15)#PCS1(15)>PCS2(15),DSURF(15)»RDQ,RDS0# F00>F01>F10,F11, 
7N0Q>lQ«KQf AWING^VOLDRG>IORGPLT( 129^15 )>SECORG( 15) 

COMMON /FLO/ STR IP^ P 1^ P 2^ P3^ B ETA# FR, IR» JR# K DG, IG# JG#KG#NS# FSWEEP 
COMMON /SWP/ DXYZ(129)#GKl(129#15)#GK2a29#15)#SX(129)#SZ(129)#SXX 
1(129)#SXZ(129)#SZZ(129)#RO(129)#R1(129)#C(129)#0(129)#G10(15)#G20( 
215)#G30(15)#G40(15)#G1( 15)#G2(15)#I1# 12# K# L# NO# LX# MX# KY# MY# Tl# AAO# 
3Q1#Q2#TYAU#S1 

COMMON /DIM/ NX1#NY1#NZ1#FDIM 

BETX-.Ol 

BETY-.15 

BETZ*.l 

BSCAL«1./(1.+FDIM) 

BSCALl-1./ (2.*(1. + FDIM) ) 

Jl-2 

IF ( FMACH.GE.l. ) Jl-3 
C( Il-1)»0. 

0(11-1) «0. 

DO 10 I»I1#I2 
R0(I)-1. 

R1(I)-1. 

GK1(I#1)-G(I#1#L) 

10 GKK I#J1-1)-6(I#J1-1#L) 

J-Jl 

I3-I2 

20 BC — Tl + BK J )*C1(K ) 

DO 60 I«I1#I3 
AB — Tl + AK I) + B1( J ) 

AC-Tl+AK I)*C1(K ) 

YP»S0(I#K)+B0{ J) 

A-1.-R0( I)+A0( I)+A0( D+YP+YP 

H»RO(I)/A 

FH«RO(I )*A 

P«AO( I )*(A,*YP*YP-FH) 

Q-YP*( A.*AO( I )*A0( I )-FH) 

A»XZ(K)tXZ(K)-YZ(K)*YZ(K) 

B-(XZ(K)+XZ(K ) )*YZ(K) 

AZ — A0( I )*XZ(K)-YP*YZ(K) 

8Z — A0(I )*YZ(K)+YP*XZ(K) 

CZ-H*H*( P+A-Q*B)“AO( I)+XZZ(K)-YP+YZZ(K) 

DZ-H+H* (Q*A+P*B)-AO( I )*YZZ(K)+YP*XZZ (K ) 

DGI»G(I+1#J#L )-G(I-l#J#L) 

DGJ»G( I#J+1#L )-GKl( I# J-1) 

DGK«G( I#J#L + 1 )-GKl(I#J) 

DGII=G( I + l# J# L)-6(I#J#L )-G(I# J#L)+6( I-1#J#L )+A3( I)*DGI 
DGJJ-G( I# J + 1# L)-G(I# J#L)-G(I# J#L)+G(I#J-1#L )-B3< J) + DGJ 
DGKK-G( I# J#L + 1)-G(I#J#L)-G(I# J# L)+G( I # J # L-1 ) +C3 ( K ) ♦DGK 
DGI J»G( I+l# J+1#L )-G(I-l# J+1#L)-G( I+l# J-1#L ) +G ( I-l# J-1# L ) 
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DGIK=G(I+1^J>L+1)-G(1+1#J>L-1 J/L+l)+G(I-i>J>L-l) 

DGJK»G( I>J+1# L + l)-G( I# J-1#L + 1 )-G ( I > J +1, L-1 ) +G ( I # J-1 , L-1 ) 

GX-AK I)*OGI 
GY—6K J )+DGJ 

U»GX-SX( I)*GY+CA*AO( I)+SA^YP 
V«GY+SA*AO{I)-CA*YP 

W*RO(I )*(C1{K)*DGK-S2(I ) *G Y+S Y AW+C A* X Z ( K ) + S A + YZ ( K ) + ( U* A Z + V*B Z ) ) 

AU-U+W+AZ 

AV-V+W+BZ 

QXY»H*(U*U+V*V) 

UU«QXY+W*W 

AA-DIM(AA0>.2+QQ) 

HZ»AZ+SX(I )“BZ+FH+SZ(n 
FXX»1.+H*AZ*AZ 
FYY-1.+SX( I)*SX( D+H+HZ+HZ 
FXY*SX( I >+H*AZ*HZ 
BV«AV-AU*SX( I )-FH*W*SZ( I) 

UU-H+AU* AU 

VV»H*BV*BV 

WW-FH+W+W 

UV-H+AU+BV 

UW>AU’i‘W 

VW-BV+W 

AXX»R1(I)*(FXX*AA-UU) 

AZZ»FH*AA-WW 

AXZ-(RO( I)+RO( I) )<‘IAZ*AA-UW) 

R— (AXX*SXX(I)+AZZ+SZZ(I)+AX2^SXZ(I))*GY+T1»(R0(I)^AA*(CZ»GX+(DZ-S 
IXd )*CZ)*GY)-H*(CA*(AU*AU-AV*AV) + (SA + SA)*AU*AV-QXY*(U*A0(I )+V*YP+( 
2W+W)*(A0(I)+AZ+YP*BZ)))-WW*(CA*XZZ(K)+SA*YZZ(K))-W+W+(U+CZ+V*DZ)) 
AXT«ABS(AU + A1 (I) ) 

AYT-ABS ( BV*B1 ( J ) ) 

AZT-ABS ( FH + W+CKK ) ) 

A«RO(I)*BETA+AA/AP,AXl(AXT»AYT>AZTf (l.-RO(l) ) ) 

AXT- A+AXT 
AYT*A*AYT 
AZT» A*AZT 

IF (UQ.GE.AA) go to 30 
AXX-AXX*A2( I) 

AYY»(FYY*AA-VV)’«'B2( J ) 

AIZ-AZZ+C2(K) 

AXY-~R1( I)+(FXY*AA+UV)*(AB+AB) 

AXZ-AXZ+AC 

AYZ«-RO( I )*(HZ*AA + VW) + ( BC + BC ) 

BP-AXX 

BM-AXX 

B — AXX“AXX“U1*(AYY + AZZ) 

R-AXX*DGI1 + AYY+DGJ J + AZZ + DGKK + AXY + DGIJ + AYZ + DGJK+AXZ + DGIK + R 
GO TO AO 
30 NS«NS+1 

S«SIGN(1.#U) 

IM»I-IFIX(S) 

IMM«IM-IFIX(S ) 
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AXX=UU*A2( I ) 

AYY-VV*B2( J ) 

AZZ«WW’4‘C2(K) 

AXY»6.*S*UV*AB 

AXZ»6.*S*UW*AC 

AYZ«8.*VW*BC 

BXX« ( FXX*QQ-UU)*A2 ( I ) 

BYY«(FYY*QQ-VV)*B2(J ) 

BZZ»(FH+QQ-WW ) + C2(K) 

BXY»-(FXY*UQ+UV)* (AB+AB) 

BXZ»( AZ*QQ-UW)*( AC+AC ) 

BYZ»-(HZ*QQ+VW)*(BC+BC) 

AQ»AA/QQ 

DELTAG»BXX*DGII+BYY*DGJ J+BZZ»D6KK+BXY*DGI J+BYZ+DGJK+BXZ+DGIK 
DGII»G(I>J>L)-G(IM>J>L)-G(IMfJj»L)+G(IMr",#J,L)+A3(I)*DGI 
OGJ J=G{ I* JjL) “Gd, J-1#L )-G( I> J-1>L) +GKK I> J-2 )-B3( J )*DGJ 
DGKK*G{ I#J,L)-G( L-l)-G( I# J# L-1)+GK2( I# J )+C3(K)*DGK 
DGIJ»G(I>J>L)-G(IM>J^L)“G( Ij J-l> L ) +G ( IK, J-1, L ) 
DGIK*=G(I,J,L)-G(I,J,L-1)-G(IK,J,L)+G(IK,J, L-1 ) 

DGJK»G( I, J, L)-G( I, J, L-1 )-G( I, J-1,,L)+G( 1, ) 

GSS«AXX*DGII+AYY»DGJ J+AZZ+0GKK4 AXY*DGIJ+AYZ*D6JK+AXZ*D6IK 
B«.5*(AU-1,)*(AXX+AXX+AXY+AXZ) 

BP-AQ+BXX-(1.-S)*B 
BM«AQ*BXX-{ 1. +S )*B 

B*-AQ*( BXX+BXX+Q2*(BYY+8ZZ) ) + ( A Q-1 . ) ♦ ( 2 . * ( A XX+ A Y Y+ A ZZ ) + AX Y + A Y Z + A XZ 
1) 

R-( AQ-1. )+GSS+AQ*DELTAG+R 
40 IF (ABS(R).LE.ABS(FR)) GQ TO i>0 
FR»R 
IR-I 
JR»J 
Kk«K 

50 R«R-AYT*(GK1( I , J - 1 ) -G { I , J -1, L) ) -AZT* ( GK 1 ( I , J )-G ( I , J , L-1 ) ) 

B-B -AXT-AYT-AZT 
BM-BK+AXT 

B«l. /{ B-BM^C ( I-l) ) 

C(I)»B*3P 

60 D( I ) »B* ( R-BM*D( I-l ) ) 

C G*0. 

I«I3 

DO 8C M»I1,I3 
CG«D(1)-C(I)*CG 

IF (ABS(CG}.L£.ABS(DG)) GO TO 70 

DG-CG 

IG-I 

JG*J 

KG-K 

70 GK2( I,J )*GK1( I,J ) 

GK1(I,J)=G{I, J,L) 

G(I,J,L)«G(I,J,L)-CG 
80 I«I-1 
J-J + 1 
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IF (J-KY) 2C>90#110 
90 IF { I2.GT. ITE2(K) ) I3*ITE2(K) 

IF ( ITE2(K) .EQ.MX) I3-LX 
DO 100 I»I1>I3 
LV=«lABS(i-IABS(IV(I^K))) 

RO(I)«AMINO(LV»IABS( IV( liK) ) ) 

100 R1(I)-LV 
GO TO 20 
110 N*NO 
I-LX+1 

IF (K.LT.KTEl .DR.K.GT.KTE2) GO TO 130 

I0-NX+2-I3 

DO 120 I>I0j13 

A«1.-R0( I)+AO(I)*AO( I)+SO{ I/K)+SO(I#K) 

H-RO( D/A 
FH»RO( I ) *A 

AZ — A0( I ) + XZ(K)-SO( I/K)+YZ (K) 

BZ»-AO(I)*YZ(K)+SO(DK)*XZ(K) 

HZ-AZ + SX(I )-BZ+FH*SZ(I ) 

FYY-l. + SX I I )*SX{ I )+H>«‘HZ+HZ 

FXY»SX( I )+H*AZ+HZ 

D6I = G( I + 1>KY, L)-G( I-1>KY>L ) 

DGK»G( DKY>L+1)-GK2 ( DKY) 

V«SA*AO( I )-CA + SO ( DK) 

U*A1(I )*DGI+CA*AO(I )+SA*SO(DK) 
W»C1CK)*DGK+SYAW+CA+XZ(K)+SA*YZ(K) 

120 G( I»KY + l»L)»6{DKY“liL)+( V*{l.-H*6Z*HZ)“UtFXY-W*HZ)/(FYY*Bl(KY) ) 
I«IO 

IF ( lO.NE. ITEKK) ) GO TO 130 
E-G( I3>KY,L)-G( I0>KY,L) 

N0*N0+1 

E0(N0)»EU(N0) +P3* ( b-EO( NO ) ) 

N-NO 

130 IF (I.LE.Il) GO TO 170 
I-I-l 
E = 0. 

IF < IV( DK) .NE.l) GO TO 160 
ZZ-Z(K)-TYAW* (XC(K)+S1*A0( I)+AC(I) ) 

140 IF ( ZZ.GE .Z0( N-1 ) ) GO TO 150 
N-N-1 
GO TO 140 

150 R»(ZZ-Z0(N-1) )/ (ZO(N)-ZO(N-l) ) 

E«R*EO(N )+(l.-R)*EO(N-l) 

160 M-NX+2-I 

G(DKY+liL)*G (M/KY-1#L )-E 
G(M>KY+l>L)-G(I^KY-l^L)+£ 

GK2(M,KY)«GK1 (M^KY) 

GK1(M>KY)-G(M>KY,L) 

G(M#KY#L)»G(I#KY#L)+£ 

GO TO 130 

C ACCORATE TRONCATED BQONDARY CONDITIONS 
170 CONTINOE 
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DO 180 I«2,NX 

180 G(I> J 1-1# U = (1.-BETY/B SCALD ♦&( I#J1#L) 
DO 190 J-l/MY 

G( Il-l# J# L ) »{ l.-BETX/BSCAL)+G( 11# J#L) 
190 G(I2+1# J#L) =( l.-BETX/BSCAL )+6(I2#J#L) 
RETURN 
END 


SUBROUTINE VELO ( K# L # S V # S M#C P #X # Y# UC # VC # WC ) 

C CALCULATES SURFACE VELUCITY 

COMMON G(12 9#12#15)#S0(12 9#i:?)#E0{131)#20(131)#IV(129#15)#ITEl(15) 
1# ITE2(15)# A0( 129)#A1{129)#A2(129)#A3{129)#B0{12)#61(12)#62(12)#B3( 
212)#Z(15)#Cl(15)#C2(15)#C3(15)#XC(15)#XZ(l!j)#XZZ(15)#YC(15)#YZ(15) 
3# YZZ( 15 ) #NX#NY#NZ#KTE1#KTE2#ISYM#KSYM#SCAL#SCALZ# YAW#CYAw# SYAw#ALP 
^HA#CA#SA#FMACH#Nl#N2#N3#I0#N0ES#TSTEP#EPSl#QPRE(129#15)#S0PRb(129# 
515)#NQSTA#2QSTA(15)#PCQ1(15)#PCQ2(15)#PCQ3(15)#QQ1( 15) ,QQ2(15)#QQ3 
6{15)#Q04(15)#PCS1(15)#PCS2(15)#DSURF(15)#RDQ#RDS0#F00#F01#F10#F11# 
7NDQ# IQ#KQ# AWING# VOLDRG#IDRGPLT( 129# 15) #SECDRG(15 ) 

DIMENSION SV(D# SM(1)# CP(D# X(D# Y(l)# UC(D# VC(D# WC(1) 
DIMENSION 02(129)# Q2S(129)# Q2SX(129)# Q2SZ(129)# Q2M(129)# Q2SM( 
1129)# Q2SXM(129)# Q2SZM(129)# Q2P(129) 

DIMENSION 0Y(129)# DYK(129) 

COMMON /DIM/ NX1#NY1#NZ1#FDIM 
I10*ITE1(KTE1 ) 

I20-1TE2(KTE1 ) 

DO 1C KDUM«KTE1#KTE2 
I10=MIN0(I1C#ITE1(KDUM)) 

10 I20«MAX0( 120# ITE2 (KDUM) ) 

J*NY+1 

Ola. 2*FMACH*+2 
Tl«l./{ .7=»FMACH* + 2) 

DO 20 I»I10#I20 

FH»AG(I)*A0(I)+S0(I#K)«S0(I#K) 

H«0. 

IF ( IV( I#K) .NE.O) H«1./FH 
AZ»-AO(I )*XZ(K)-SO(I#K)*YZ{K) 

BZ«-AO( I)*YZ(K)+SO( I#K)*XZ(K) 

DSI-SO( I+1#K)-S0( I-1#K) 

DSK»S0(I#K+1)-S0(I#K-1) 

SX«A1(I)*DSI 
SZ«C1(K )^DSK 

DGI«G(I + 1# J#L )-G( I-1#J#L ) 

DGJ«G(I#J+1#L)-G(I#J-1#L) 

DGK*G{ I# J# L + 1 )-G( I#J#L-1) 

U»A1(I)*DG1+SX*B1( J )*DGJ+CA*A0( I )+SA*SC(I#K) 
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J )*DGJ+SA+AO< I )-CA*SO( I#K ) 

W-CKKj + DGK+SZ+BKJ )*0GJ + SYAW+CA*XZ(K) + SA+YZ(K)+H+(U*AZ + V*B2) 

QQ«H+(U+U+V*V)+W*W 

SV( I ) =SIGN<SQRT{QQ) #U) 

IF ( IV( I^K) .EQ.O) SV(I)-SViI-l)+SV(I-l)-SV(I-2) 

UC( I)-0. 

VC (I )»0. 

WC(I)-0. 

QQ-1.+Q1*(1.-QQ) 

SM(I ) «FMACH+SV(I )/SQRT ( QQ) 

CP( I)»Tl+(QQ+*3.5-l. ) 

X(I)«XC(K)+.5^SCAL*( AO( I)*AO( I)“SO(I>K)+SO{IjK) ) 
Y(I)«YC(K)+SCAL*AO(I )*SO(I>K) 

IF (NDES.LE.O) GO TO 20 

W2S-2.*W*H*( {-U*YZ(K)+V*XZ(K)+SA*AZ-CA*BZ)-2.*S0( I>K)*H*(U»AZ+V*BZ 
1 ) ) 

Q2S ( I ( (SA + U-CA + V )-S0( I#K)*H» (U+U + V + V ) )-w2S 
Q2fX ( I ) =“2.*B1( J )*OGJ*H*(U+AZ*W) 

Q2SZ( I) »-2.*Bl( J)*OGJ*W 
20 CONTINUE 

IF (NDES.LE.O) RETURN 

CALL SURMOO ( Kf L^ SV# SM> CP# X, Y^ UC# VC# WC ^ Q2S> Q2SX> Q2SZ ) 

RETURN 

END 


SUBROUTINE SURMGD ( K # L # S V # SM , C P # X# Y # UC # VC # l«C # Q2 S # 0 2S X # Q2S Z ) 

C PERFORMS SURFACE MODIFICATION IN DESIGN MODE 

COMMON G(129#12#15)#S0(129#15)#E0(i31)#Z0{131)#IV(129#15)#ITEl(15) 
l#ITE2(15)#AO( 129)#Al{129)#A2(129)#A3(lZ9)#B0(12)#Bi(12)#ti2(12)#B3( 
212)#Z(15)#C1(15)#C2(15)#C3(15)#XC(15)#XZ(15)#XZZ(15)#YC(15)#YZ(15) 
3#YZZ{15)#NX#NY#NZ#KTE1#KTE2#ISYM#KSYM»SCAL#SCALZ#YAW#CYAW#SYAW#ALP 
^HA#CA#SA#FMACH#Nl#N2#N3#ia#NDES#TSTEP#EPSl#QPRE( 129#15)# S0PRE( 129# 
51S)#NQSTA#ZQSTA(15)#PCQl(15)#PCQ2(15)#PCQ3(lb)#QQl(15)#QQ2(lD)#QQ3 
6{15)#QQ4(15)#PCS1(15)#PCS2(15)#DSURF(15)#RDQ#RDS0#F00#F01#F10#F11# 
7NDQ#IQ#KU# AWING# VOLDRG#IDRGPLT( 129# 15)#S£CDRG(15) 

DIMENSION SV(D# SM(D# CP(D# X(D# Y(D# UC(D# VC(D# WC(1) 
DIMENSION 02(129)# Q2S(129)# Q2SX(129)# Q2SZ(129)# Q2M(129)# 02SM( 
1129)# Q2SXM(129)# Q2SZM(129)# Q2P(129) 

DIMENSION DY(129)# DYK(129) 

COMMON /DIM/ NX1#NY1#NZ1# FDIM 
I1=ITE1(K) 

I2»ITE2(K) 

IlO-ITEKKTEl ) 

I20«ITE2(KTE1) 

DO 10 KDUM«KTE1#KTE2 
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I10*MIN0{I10j IT fcl (KDUM) ) 

10 I20-MAX0{ I20> ITE2(KD0M) ) 

DX»1. /FLOAT(NXl) 

CUTL0«.80 

J*NY+1 

Ql».2*FMACH+*2 

Ti-l./(.7*FMACH*+2) 

KFLAG-0 

IF (K.GT.KTEl) GO TO 30 
DU 20 I=I10#I20 
DYK( I) -0. 

Q2(I)«QPRE(I#K)*QPRE(1>K)-SV(I)=4‘SV(I) 

UC (I ) *QPRb ( I>K) 

VC(I)-SO(I#K)-SOPRE(I>K) 

Q2SM( I)»Q2S(I ) 

Q2SXM( I ) =02SX{ I ) 

20 Q2SZM ( I ) -Q2SZ ( I ) 

RETURN 

30 DO ^0 I»I10,I20 

Q2P(I)='0PRE{l>K)*gPRE(i>K)-SV{I)=«^5V(I) 

UC(I)»QPRE(I#K) 

40 VC ( I ) »S0( I>K )-SOPRE { I,K ) 

IF (K.GT.KTEl+1) GO TO 60 
DO 50 I=I10>I20 
50 Q2M( I ) »Q2P (I ) 

60 CONTINUE 
DT»TSTEP*DX 
KM=>K-1 

IF (KFLAG.EQ.l) KM«=KT£2 

I«I2C 

DY(I ) =0. 

DYK( I) -0. 

70 I-I-l 
DY(I )=0. 

IF ( I.LE. 110) GO TO 100 

IP-I+1 

IM*I-1 

IF (I .GT.ITE2(KM) .QR.I.LT.ITEKKM) ) GO TO 70 
Q2X*2.*Ai(I)+(Q2(IP) -0 2(D) 

IF ( Q2SXM( I ) . LE.O. ) 02 X -2 . ♦ A1 ( I ) + ( 02 ( I ) -02 ( IM ) ) 
Q2Z«2.*C1(KM)*(Q2P(I)-02(I)) 

IF (02SZM( I) . LE.O. ) 0 2 Z -2 . *C 1 ( K M ) * ( 0 2 ( I ) -0 2M ( I ) ) 

DS»02(I ) + .5*DT*(02SM(I )«Q2(D+Q2SXM( D*Q2X+02SZM( I )*02Z) 

FFOO-FOO 

FFOl-FOl 

FFIO'FIO 

IF { ABS(SN ( I) ) .LE.CUTLQ) FFlO-0. 

FFll-Fll 

FaC«1./(FF00+FF0I-FF10-FF11) 

DY(I ) «( DT*DS-OY( IP)+(FF10+FF11)+DYK( I )♦ (FF01-FF11)+DYK( IP)*FF11)*F 
lAC 

DUMM»SO( DKH) +DY( I ) 
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IF (DUMM.LT.SOPRE{I#KM) ) DUMM »S OPRE ( I > K M ) 
DY( I) -DUMM-S0( I^KM) 

SO( I>KM ) -DUHM 

IF (SO(I#KM).LE.SOPRE(I#KM) ) GO TO 90 
IF ( ABS( DS ) .LT.ROSO) GO TO 80 
RDSO-AMAXKRDSO^ ABS{DS) ) 

IQ«I 

KQ«KM 

80 CONTINUE 
NDQ«NOQ-i-l 
RDQ»RDQ+A6S(DS) 

90 CONTINUE 
GO TO 70 
100 CONTINUE 

IF (KFLAG.EQ.l) RETURN 
DO 110 1-110,120 
DYK(I)-DY(I) 

Q2M(I)-Q2(I) 

Q2( I ) -Q2P ( I ) 
a2SM(I)»Q2S(I ) 

Q2SXM( I ) =U2SX { I ) 

Q2SZM( I ) -U2SZ ( I ) 

110 CONTINUE 

IF (K.LT.KFEZ) RETURN 
DO 120 1-110,120 

120 Q2P( I)»AMIN1(0.,2.*Q2(I)-Q2M(I) ) 

KFLAG-1 
GO TO 60 
END 


SUBROUTINE DRAGC (IND,SCALX) 

C COMPUTES THE WAVE DRAG BY VOLUME INTEGRATION OF ENTROPY INEQUALITY 
COMMON G( 12 9, 12, 1 5 ) , S 0 ( 129, 15 ) , E 0 ( 1 31 ) , ZO ( 1 3 1 ) , I V ( 1 29 , 15),ITE1(15) 
1,ITE2(15),A0(129),A1(129),A2(129),A3(129),130(12),B1(12),B2(12),B3( 
2i2),Z(15),Cl(15),C2(15),C3(15),XC(i5),XZ(i5),XZZ(15),YC(15),YZ(15) 
3,YZZ(15),NX,NY,NZ,KTE1,KTE2,ISYM,KSYM,SCAL,SCALZ#YAW,CYAW,SYAW,ALP 
^HA,CA,SA,FMACH,N1,N2,N3,I0,NDES, TSTEP,EPS1,QPRE ( 129, 1 5 ) , SOPR E ( 1 29, 
515),NQSTA,ZQSTA( 15),PCQ1(15),PCQ2(15),PCQ3(15),QQ1(15),QQ2(15),QQ3 
6(15), QQ ^(15), PCS1(15),PCS2(15),DSURF( 15),kDQ,RDS0,F00,F01,F10,Fll, 
7NDQ, IQ,KQ, AWING,VOLDRG, IDRGPLK 129,15),SEC0RG(15) 

COMMON /FLO/ 5TRIP,Pl,P2,P3,BETA,FR,IR,Jk,KR,DG,lG,JG,KG,N5,FSWEEP 
COMMON /SWP/ DXYZ(129),6K1(129, 15),GK2(129,15),SX(129),SZ(129),SXX 
1(129),SX2(129),SZ2(129),R0(129),R1(129),C(129),D(129),G10(15),G20( 
215),G30(i5),G40(15),Gl(15),G2(15),Il,I2,K,L,N0,LX,MX,KY,MY,Tl, AAO, 
3Q1,Q2,TYAW,S1 
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DIMENSION XXX(3^), YYY(3<r), DRAG(3^> 

COMMON /DIM/ NX1,NY1,N21>FDIM 

IF { SCALX.EQ.0.0) SC AL X»5 . / ( Z (KTE2 )-Z (KTE 1 ) ) 

TTTX-3. 

SSSX«-SCALX*XC(KTE1 ) 

DX=2./FL0AT(NX1) 

DZO-l./FLOAT<NZl) 

DVOL-2./ FLOAT (NX1*N Yl+NZl ) 

DO 10 K»KTE1#KTE2 
10 SECDRG(K)«0. 

FACZ-.5 

LX«NX/2+l 

PI-3.1415927 

RAO-PI/180. 

ANG — 5.*RAD 
NDARK-40 
SCOT-. 02 
INDPLT-0 

SZ0-200.+ (FLOAT(NX) /156. )**2 
DO 130 K»KTE1>KTE2 
INDPLT-INDPLT+1 
DO 20 J-2#KY 
XXX( J )«0. 

YYY( J ) *0. 

20 DRAG(J)«0. 

SSSY»5.*(Z(K)-Z(KTE1) )/(Z(KTE2)-Z(KTEl) )+2.45 

KP-K+1 

KM-K-1 

I1«ITE1(K) 

I2-ITE2(K) 

DO 50 I»I1#I2 
IDRGPLK I>K)«-1 
IP-I+1 
IM-I-1 

SX( I )«A1 ( I )♦( S0( IP# K)-S0( IM#K ) ) 
SZ(I)*C1(K)*(SO(I#KP)-SO(I#KM)) 

FACY-1. 

DO 40 J«2#KY 

IF (J.EQ.KY) FACY-.5 

YP«SO(I#K)+BO(J) 

FH-A0( I)+AO(I )+YP+YP 
H-l./FH 

AZ — A0( I ) + XZ<K)-YP+YZ(K) 

BZ — A0( I)*YZ(K)+YP*XZ(K) 

DGI»G(IP#J#K)-G(IM#J#K) 

DGJ«G( I# J+1#K )-G( I# J-1#K) 

DGK«G(I#J#KP)-G(I#J#KM) 

GX-AKI )*DGI 
GY--BK J ) + OGJ 

U-GX-SX( I )*GY+CA+A0( 1)+SA*YP 
V-GY+SA*AO(I)-CA*YP 

W-(C1(K ) + DGK-SZ( I )*GY + SYAW+CA*XZ(K) +SA»YZ(K)+H=«' (U+AZ + V + BZ) ) 
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AU-U+W+AZ 

QXY«H*{U+U+V*V) 

QQ-QXY+W+W 
AA-01M( AAOf *Z*QQ) 

DUMM«0. 

IF (OQ.LT.AA) GO TO 30 
UU-H+AU+AU 
AXX«UU*A2( I) 

AQ-AA/QQ 

XJAC0«1.*SCALZ*( ( l.+4.*B0( J)*BO( J) )**1.5)*DVOL/SCAL 
DRGSS»(G( I+l> J^K)~2.*G( J,K)+G( I-l# J,K ) ) / (DX+DX) 
DRGSS-.5*(0RGSS-ABS(DRGSS) ) 

ROCS*FMACH+FMACH* ( FM ACH*F M ACH*A A ) ♦♦! . 5 

OUI1M»2.»DX#SC AL*SCAL+( l.-AQ)*RDCS*AXX*DRGSS*ORGSS+SQRT (H )^XJACD/AW 
lING 

C SECOND ORDER ACCURATE VOLUME INTEGRAL 
VOLDR6»VOLORG+FACY*FACZ*OUMM 
SECDRGCK ) «SECDRG(K) +AWING+DUMM/ (SCALZ^DZO) 

IF (J.EQ.KY) IDRGPLT( I#K)-1000000.*DUMM 
30 IF (IND.NE.l) GO TO ^0 
IF (I.LT.LX) GO TO 40 

XX-XC(K) + .5*SCAL+( A0( I)*A0( I)-(SO(I^K) + BO( J ) )*»2 ) 

YY-YC(K) + SCAL*AO( I ) + (S0( I>K) + BO( J) ) 

DRAG( J)«DRAG( J)+DUMM 
XXX( J)-XXX< J)+DUMM+XX 
YYY(J)»YYY(J)+OUMM*YY 
40 CONTINUE 
50 CONTINUE 

IF (IND.NE.l) GO TO 130 
DO 60 J»2^KY 

IF (DRAG( J ) .LT.l.E-50) GO TO 60 
XXX( J)«XXX( J) /DRAG{ J ) 

YYY( J)«YYY( J) /DRAG( J ) 

60 CONTINUE 
IPREV-0 

DO 120 JJ-2#KY 
J-KY+2-JJ 
SIZE«SZ0*DRAG( J ) 

IF (SIZE. LT. SCUT) SIZE-0. 

XC0«SCALX*XXX(J )+SSSX+TTTX 
YC0»SCALX*YYY( J )+SSSY 
XL-XCD-SIZE*C0S(ANG) 

YL-YCD-SI2E*SIN(ANG) 

IF (SIZE. LT. SCUT. AND. IPREV.EQ.O) GO TO 110 
IF (J.EQ.KY) GO TO 110 

IF ( IPREV.EQ.l. AND. SIZE. GE. SCUT ) GO TO 80 
IF ( SIZE. LT. SCUT) 60 TO 70 
XX-2<:*( XXX(J)-XC(K) )/SCAL 
YY-2.*(YYY( J)-YC (K) )/SCAL 
RR-SQRT(SQRT( XX*XX+YY*YY) ) 

THET-.5*ATAN( YY/XX ) 

XX-RR»COS(THET) 
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YY«RP*SIN(ThET) 

YY»YY-BO( J )+B0( J+1) 

XCDO«XC (K)+.5*SCAL*(XX*XX-YY*YY) 

YCDO«YC(K)+SCAL*XX*YY 

XCOO-SCALX+XCDO+SSSX+TTTX 

YCDO-SCALX+YCDO+SSSY 

XLO»XCiJO 

YLO-YCDO 

GO TO 80 

70 XX»2.+(XXX( J+1)-XC(K) )/SCAL 
YY-2.*(YYY< J+1)-YC(K) )/SCAL 
RR»SQRT(SQRT(XX*XX+YY*YY)) 
THET-.5+ATAN( YY/XX) 

XX-RR+COS (THET) 

YY«RR*SIN(THET) 

YY«YY-BO( J+1)+B0(J ) 

XCO»XC (K)+.5*SCAL*(XX*XX-YY*YY) 
YCD«YC(K)+SCAL»XX*YY 
XCO-SCALX+XCD+SSSX+TTTX 
YCD»SCALX*YCD+SSSY 
XL-XCD 
YL-YCD 
80 CONTINUE 

IF (KTE2.lt. 10) GO TO 90 
IF (MOD( INDPLT^2) .EQ.O) GO TO 110 
90 CONTINUE 

DO 100 L*1jNDARK 
FAC»FLOAT(L )/ FLOAT (NO ARK ) 
XX«FAC*XCD+(1 .-FAC)»«XCD0 
YY-FAC*YCD+(1.-FAC)+YCD0 
CALL PLOT (XX#YY#3) 
XX»FAC*XL+(1.-FAC )*XLO 
YY«FAC*YL+(1.-FAC )*YLO 
100 CALL PLOT (XXfYY#2) 

110 IPREV-0 

IF (SIZE. GE. SCUT) IPREV=1 

XCOO=XCO 

YCDO-YCD 

XLO-XL 

YLO»YL 

120 CONTINUE 
130 CONTINUE 
FACZ-1. 

RETURN 

END 
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SUBROUTINE CPLDT ( I 1# 12^ X# 2^ A# 0> FMACH ) 

C PLOTS CP AT EQUAL INTERVALS IN THE MAPPED PLANE 

DIMENSION K0D£(3)# LINE(IOO)# X(l)» Yd), A(l), B(l), C(l)> D(l)> 
1Z(1) 

DATA KOOE/IH ,1H+,1H0/ 

IWRIT=6 

AA0«.2+1./FMACH**2 
WRITE (IWRIT#80) 

DO 10 1-1,100 
10 LINE(I)-KUDE(1) 
liO-Il 
120-12 
ISPA-1 

LX-( 11+12)22 
FDEN-1. /B( 120 ) 

LX0».85*FL0AT(LX) 

IF < LXO.LT. 12-^9) ISPA-2 
AMAX-0. 

CMAX-0. 

DO 20 1-110,120 
AMAX»AMAX1(AMAX,ABS(X(I))) 

AMAX-AMAX1(AMAX,ABS(Y(D) ) 

CMAX-AMAX1(CMAX,ABS(C( I) ) ) 

20 CMAX-AMAXKCM AX, ABS (D( I ) ) ) 

DO 70 I»LX0,I2,ISPA 
XFRAC»FDEN + B( I ) 

Y(I)-SQRT(Y(I )**2/(AA0-.2+Y(I)**2) ) 

Kl-(59./AMAX)*ABS(X(I))+4il. 

K1-MIN0(K1,100) 

LINE (Ki ) -KODE (2) 

K2-(59./AMAX)*ABS(Y(I))+41. 

K2-MIN0(K2,100) 

LINE(K2)«KUDE (3) 

K3=( 36. /CMAX) *0( I )+l. 

K3-MIN0(K3,40 ) 

IF (K3.GE.1) GO TO 30 
K3-1. 

GO TO AO 

30 LINE (K3) =KODE { 2) 

AO KA-(36./CMAX)#C(I)+l. 

IF (KA.GE.l) GO TO 50 
KA-1 

GO TO 60 

50 LINE(KA)»KODE (3) 

60 CONTINUE 
J J-0 

IF (Z(I).GT.O) JJ-1 

WRITE (IWRIT,90) XF R AC , X ( I ) , Y ( I ) , J J , L IN E 
LINE(K1 ) -KODE (1) 

LINE(K2)-K0DE(1) 

LIN6(K3)-K0DE (1) 

LINE(KA) -KODE(l) 


69 



70 CONTINUE 

11- IlO 

12 - 120 
RETURN 

80 FORMAT ( IX, /, 3X> 3HX/C > 5HMC0MP# AX# 5HMDSGN # 3H ON) 
90 FORMAT ( IX, F5 . 2, 2F9. 3, I 3, 2X, lOOAl) 

END 


SUBROUTINE FORCE ( 1 1, 1 2, X, Y,C P, AL,C HORD, XM, CL ,C D, CM ) 
C CALCULATES SECTION FORCE COEFFICIENTS 
DIMENSION X(l), Y(l), CP(1) 

RAD- 57. 2957795130823 

ALPHA-AL/RAO 

CL-O. 

CD-0. 

CM-0. 

N-I2-1 

DO 10 I»I1,N 
DX-( X( I+l )-X( I ) ) /CHORD 
DY-(Y( I+l)-Y( I) )/CHORD 
XA-(.5*{X(I+1)+X(I) )-XM) /CHORD 
YA-, 5*{ Y( I+l) +Y ( I ) ) /CHORD 
CPA».5+(CP (I + 1)+CP(I ) ) 

DCL— CPA*DX 
DCD-CP A+DY 
CL-CL+DCL 
CD-CD+DCD 

10 CM»CM+DCD*YA-DCL*XA 

DCL-CL+COS (ALPHA ) -CD* SIN (ALPHA) 
CD-CL*SIN(ALPHA)+CD*CQS( ALPHA) 

CL-DCL 

RETURN 

END 


SUBROUT IN E TOT FOR ( KTE 1,K TE2, CHORD, SCL,SCD, SCM, Z, XC,CL, CD, CMP, CMR, 
ICMY, AWING) 

C CALCULATES TOTAL FORCE COEFFICIENTS 
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DIMENSION CHGRD(1)> SCL(1)> SCD(l)/ SCM(1)> Z(l)> XC(1) 
SPAN»Z(KTE2)-Z(KTE1) 

CL-0. 

CO-0. 

CMP-0. 

CMR-0. 

CMY«0. 

S»0. 

N-KTE2-1 
DO 10 K-KTEl^N 
DZ-.5*{Z(K + 1)-Z(K) ) 

AZ-.5+(Z(K+l)+Z(K)) 

C L-C L + DZ*(SCL(K + 1)*CH0RD(K + 1)+SCL(K) + CHORD (K) ) 

CO«CD+DZ+(SCD (K+1)*CHQRD(K+1) +SCD(K)+ChQRD (K ) ) 

CMP-CMP+DZ* (CHORD (K+D+ (SCM(K+1 )+CHQRD(K+l )-SCL (K+1 )tXC (K+1) )+CHOR 
1D(K)*(SCM(K)*CH0RD(K)-SCL(K)*XC (K) ) ) 

CMR-CMR+AZ*DZ*( SCL (K+1 )*CH0RD<K+1)+SCL(K)*CHQRD (K) ) 

CMY«CMY + A2+DZ>»‘( SCD(K+1)*CH0RD(K+1>+SCD(K)*CH0RD(K> ) 

10 S«S+D2*{CHQRD(K+1)+CHQRD(K) ) 

AWING«S 

CL-CL/S 

CD-CO/S 

CMP-CMP+SP AN/ S**2 
CMR»(CMR+CMR) /(S*SPAN) 

CMY»(CMY+CMY) /(S*SPAN) 

RETURN 

END 


SUBROUTINE OPT ( UQl^ 0U2> QU3# QQ^ ) 

C INITIALIZES PARAMETERS FOR THE OPTIMIZATION ROUTINE 
C AND CALLS OPTIMIZER 

DIMENSION QQ1(1)> QQ2(l)i DG3(1)# QQA(l) 

DIMENSION X(ZO), G(ZO)f H(20»20)> W(bO)> XM(20) 

EXTERNAL ORFCT 

IWRIT-6 

N»4 

DFN«.0003 

HH«1, 

MODE-1 
DO 10 I-1>N 
10 XM(I)-.05 
X(1)»QQ2(2) 

X(2)*QQ3(2) 

X(3) -QQ2(3) 

X ( 4) -QQ3 { 3 ) 
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MAXFN*30 
IPRIM = 1 
EPS-.l 

CALL VAIOA (DRFCT^N#X^ H#W/DFN» XM#HH#EPS/MODE>MAXFN^ IPRINT, lEXI 
IT) 

RETURN 

END 


SUBROUTINE DRFCT (NDUM> XOUM, F ) 

C EVALUATES THE DRAG AS A FUNCTION OF THE SPEED DISTRIBUTION 
C DETERMINED BY THE OPTIMIZER 

COMMON G(129, 12>15)>S0(129j15)>E0(131)>Z0(131)jIV{ 129,1B),ITE1(15) 
1>ITE2(15), A0(129)»A1(129),A2{129)#A3(129),B0(12) >Bl(12)>B2a2)>B3( 
212),Z(15)^Cl(i5)#C2(15),C3{15),XC(i5),XZ(15)#XZZ(15)#YC(i5),YZ(15) 
3> YZZ( 15)j»NX,NY^NZ>KT61>KTE2>ISYM>KSYMjSCAL#SCALZ#YAWf CYAW^SYAW^ ALP 
AHA,CA, SA,FMACH,Nl,N2/N3,IO»NDES>TSTEP,EPSl>QPRE(129,15)>SOPRE(129, 
51f>)>NQSTA/ZQSTA(15)> PCQl(15),PCQ.2(15)>PCa3( 15)>QU1( lb)>QQ2(15)#QQ3 
6(15)^QQ4a5)>PCSl(15>),PCS2(15),DSURF(15)#RDQ^ RDSO,FOOf F01,F10,F11, 
7NDQ>1U>KQ#AWING,V0L0 RGjIDRGPLT( 129, 15 ) , SEC DR6 ( 1 5 ) 

COMMON /FLO/ STRIP,P1,P2,P3,B6TA,FR,IR,JR,KR,DG,IG,JG,KG,NS,FSWEEP 
DIMENSION XDUM(l) 

dimension SV(129), SM(129), CP(129), X(129), Y(129), UC(129), VC(1 
129), WC(129), CHORDdS), SCL(15), SCD(15), SCM(15) 

DATA ISTART/l/,NlTQT/0/ 

DATA VAR/0.0/ 

IWRIT-6 

AA0»1./FMACH**2+.2 

NE-129 

LX-NX/2+1 

QCUT- .015 

UMCUT».09 

IMAX«150 

QQ2(2)»XDUM(1 ) 

UQ3(2)«XDUM(2) 

QQ2(3)«X0UM(3) 

QQ3(3)»XDUM{^) 

IF (ISTART.NE.l) CALL TREAD 
ISTART»ISTART+1 
IF ( ISTART.GE .10) IMAX-200 
DO 10 KL-1,5 

PRINT 70, ZQSTA(KL),QQl(KL),aQ2(KL),0Q3(KL),QQ4(KL) 

10 CONTINUE 

CALL SETQS (NE,NX,QPRE,S0,S0PRfc,ITEl,ITE2,KTEl,KTF2,Z,ZQSTA,A0,PCQ 
li,PC02,PCQ3,UC,VC,Q01,UQ2,QQ3,QQ^,PCSl,PCS2,DSURF,NQSTA) 

WRITE (IWRIT,120) 
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WRITE (IWRIT,80) 

ITER«0 

20 ITER»ITER+1 
CALL MIXFLO 
NIT0T«NIT0T+1 
VOLDRG-0. 

CALL DRAGC (0) 

RDS0»0. 

NDQ«0 

RDQ»0. 

DO 30 K=KTE1,KTE2 

CALL VELQ ( K, S V# SM, C Y# LC> VC# WC ) 

Il-ITEKK) 

I2»ITE2(K) 

CHORD(K)«X( Il)-X(LX) 

CALL FORCE (I1#I2#X#Y#CP#AL#CH0RD(K)#XC(K)#SCL(K>#SCD(K)#SCM(K)) 

30 CONTINUE 

CALL TOT FOR (KTE1#KTE2# CHORD# SCL»SCD#SCM#Z#XC#CL#CD1#CMP#CMR#CMY#A 
1WIN6) 

IF (NDQ.GT.O) ROQ*RDQ/FLQAT(NDO> 

WRITE (IWRIT#90) IT E R# F R# DG# N S# RDQ# R D SO # VO L DRG# C L # N I T OT 
IF ( ITER.GE. IMAX) GO TO -^0 
IF ( AMAXK RDQ#RDSO) .GE.2. ) GO TO 40 
IF (RDQ.GT.QCUT.OR.RDSO.GT.OMCUT) GO TC 20 
40 VQLDRG-0. 

CALL DRAGC (0) 

WRITE (IWRIT#100) VOLDRG 

NDUM0*NDUM/2 

DO 50 I«2#3 

XMA1 = SQRT (QQ2 ( I )**2/ ( AAO- . 2 + QQ2 ( I ) **2 ) ) 

XMA2-SQRT(QQ3 ( I ) **2 / ( AAO- . 2*003 ( I )**2 ) ) 

50 WRITE (IWRIT#110) I # 0Q2 ( I ) # QQ 3 ( I ) # XM A 1 # XM A2 
DO 60 K»KTE1#KTE2 

CALL VELO (K#K#SV#SM#CP#X#Y#UC#VC#WC) 

11- ITEKK ) 

12- ITE2(K) 

CHORD(K) »X( Il)-X(LX) 

CALL FORCF (I1#I2#X#Y#CP#AL#CHORO(K)#XC(K)#SCL{K)#SCD(K)»SCM(K)) 
WRITE (IWRIT#120) Z(K)#SCL(K) 

CALL CPLOT (I1#I2#SM#UC#VC#0PRE(1#K)#AC#S0PRE(1#K)#S0(1#K)#FMACH) 
60 CONTINUE 
F-VOLDRG 

CALL CHEKPTX (VAR) 

RETURN 

C 

70 FORMAT (1X#5F10.5) 

80 FORMAT (1X#6X#4HITER#5X#5HR£SID#6X#4HDPHI#8X# 2HNS# 5X# 5HDQAVE# 5X# 5H 
1DQMAX#6X#4HDRAG#8X# 2HC L# 5 X# 5HN I TOT# / / / ) 

90 FORMAT ( IX # II 0# 2 E 10 . 3# I 10# 2E 10. 3# FIO . 5# F6 . 2 # I 6 ) 

100 FORMAT ( IX# // #1X5HDRAG-# F10.5#20H SPEEDl SPEED2#20H MACHl 
1 MACH2 ) 

110 FORMAT ( IX# /# 6X# 110# 4F10. 5) 
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120 FORMAT ( IX, 7HZ ( K ) = ,F10.5,2X,bH CL * ,F10.2//) 
END 


SUBROUTINE TREAD 

C READS THE POTENTIAL STORED AT THE END LF THE LAST LINE SEARCH 

COMMON G( 129, 12,15), S0( 129, 15 ), E0( 131 ) , Z0( 131),IV(129,15),ITE1(15) 
1,ITE2(15),A0(129),A1(129),A2(129),A3(129),B0(12),B1(12),B2(12),B3( 
212),Z(15),C1(15),C2(15),C3(15),XC(15),XZ{15),XZZ(15),YC(15),YZ(15) 
3,YZZ{15),NX,NY,NZ,KTE1,KTE2,ISYM,KSYM,SCAL,SCALZ,YAW,CYAW,SYAW,ALP 
AHA,CA,SA,FMACH,M,N2,N3,IQ,NDES,TSTEP,EPS1,QPRE ( 129, 15), S OPR E( 129, 
515),NQSTA,ZQSTA(15),PCQ1(15),PC02(15),PCQ3(15),QQ1(15),QQ2(15),QQ3 
6(15),QQ4{15),PCS1(15),PCS2(15),DSURF{15),RDQ,RDSO,FOO,F01,F1G,F11, 
7NDQ, IQ,KO, AWING, VOL DRG, IDRGPLT ( 129, 15 ) , SECDRG ( 15 ) 

NT«8 

REWIND NT 

READ (NT) NX,NY,NZ,NM,K1,K2,NIT 

MX-NX+1 

MY-NY+2 

MZ-NZ+3 

DO 10 K-1,MZ 

READ (NT) ( (G ( I, J,K ) , I»1,MX ), J*1,MY) 

10 CONTINUE 

READ (NT) (E0(K),K=K1,K2) 

READ (NT) ((S0(I,K),I«1,MX),K«1,MZ) 

REWIND NT 

RETURN 

END 


SUBROUTINE TWRIT 

C STORES THE POTENTIAL AT COMPLETION OF LINE SEARCH 

COMMON G(129,12,15),S0(129,15),E0(131),Z0(131),IV(129,15),ITE1(15) 
1,ITE2(15),A0( 129),A1(129),A2(129),A3(129),B0( 12),B1(12),82(12),B3( 
212),Z(15),C1(15),C2(15),C3(15),XC(15),XZ(15),XZZ(15),YC(15),YZ(15) 
3,YZZ(15),NX,NY,NZ,KTE1,KT£2,ISYM,KSYM,SCAL,SCALZ,YAW,CYAW,SYAW,ALP 
AHA,CA,SA,FMACH,N1,N2,N3,I0,ND£S,TSTEP,EPS1,QPRE(129,15),S0PR£(129, 
515),NQSTA,ZQSTA( 15),PC01(15),PCQ2(15),PCQ3(15),QQ1( 15),QQ2(15),0Q3 
6(15),QQ^(15),PCS1(15),PCS2(15),DSURF(15),RDQ,RDS0,F00,F01,F10,F11, 
7NDQ, IQ,KQ, AWING, VOL DRG, IDRGPLT (129, 15 ) , SEC DRG ( 1 5 ) 
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NT-8 

REWIND NT 
NM-0 
NIT-0 
Kl-KTEl-1 

K2»KTE2+ITE2(KTE2 )-NX/2 

WRITE (NT) NX,NY/NZ,NM/K1>K2^ NIT 

MX-NX+1 

MY-NY+2 

MZ-NZ+3 

DQ 10 K-1>MZ 

WRITE (NT) ((G(I>J#K),I-1,MX),J»1>MY) 
10 CONTINUE 

WRITE (NT) (E0(K),K-K1,K2) 

WRITE (NT) (( S0(IjK)>I=1>MX),K»1#MZ) 
REWIND NT 
PRINT 20 
RETURN 

20 FORMAT (1X>15H WRITE ON TAPES) 

END 


SUBROUTINE VAIOA ( FUNCT»N^X>F#6>H»W>DFN#XM>HH» EPS# MODE# MAX FN^IPRIN 
IT# lEXIT) 

C OPTIMIZATION SUBROUTINE 
C PERFORMS LINE SEARCH FOR DRAG MINIMUM 
REAL X( 1 ) # G( 1 )#H( 1) # W ( 1 )# XM (1 ) 

IF ( IPRINT.NE ,0) PRINT 190 

IF ( IPRINT.NE .0) PRINT 200# DFN#HH# ( XM ( I ) # I -1# N ) 

NN»N*(N+1) /2 

IG-N 

IGG-N+N 

IS-IGG 

IOIFF-1 

IEXIT-0 

IR-N 

IF (M0DE.EQ.3) GO TO 'tO 
IF (M0DE.EQ.2) GO TO 30 
IJ-NN+1 
DO 20 I-1#N 
DO 10 J-1#I 
IJ-IJ-1 
10 H( IJ )-0. 

20 H(IJ)«1. 

GO TO AO 
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30 CONTINUE 

CALL MCllB (H#N,IR) 

IF (IR.LT.N) RETURN 
AO CONTINUE 
Z«F 
ITN»0 

CALL FUNCT (N>X,F) 

CALL TWRIT 

IFN»1 

DF«DFN 

IF (DFN.EQ.O.) DF«F-Z 
IF (DFN.LT.O.) DF«ABS(DF*F) 

IF (DF.LE.O.) DF=1. 

50 CONTINUE 
GO TO 170 
t)0 CONTINUE 

IF ( IFN .G£ .MAXFN ) GO TO 130 

IF (IPRINT.EQ.C) GO TO 70 

IF (MDD{ ITN/ IPRINT) .NE.O) GO TO 70 

PRINT 210> ITN,IFN 

PRINT 220# F 

IF ( IPRINT. LT .0) GO TO 70 
PRINT 220# (X(I)#I«1#N) 

PRINT 220# (W(IG + D# 1 = 1#N) 

70 CONTINUE 
ITN-ITN+1 
DO 80 I«1#N 
80 W (I ) =-W ( IG+I ) 

CALL MCllE (H#N#W#G#1R) 

Z*0. 

GSO-0. 

DO yo I»1#N 

w(is+i)*w(i) 

IF (Z*XM(I).GE,ABS(W(I) ) ) GO TO 90 
Z*ABS(W( I) )/XM(I) 

90 GSO»GSO + W( IG+I )*W ( I ) 

IEXIT-2 

IF (GSO.GE.O.) GO TO lAO 
ALPHA — 2. + DF/GSC 
PRINT 230# ALPHA 
IF ( ALPHA, GT. 1. ) ALPHA-1. 
DSTEP0-,C3/SQRT(ABS (GSO) ) 

JMIN-1 
FDMIN-F 
DO 110 JD-2#5 
DSTEP»FLQAT( JD-1)*DSTEP0 
DO ICO I»1#N 

100 »J(I)-X{ I)+DSTEP*W(IS + I) 

CALL FUNCT (N#W#F1) 

IF (Fl.GE.FDMIN) GO TO 110 

FOMIN-Fl 

JMIN-JD 
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110 CONTINUE 

DSTEP»FLOAT( J MIN-1 )*DSTEPO 
DO 120 I=1>N 

120 W(I)»X( I)+DSTEP+W(IS+I) 

CALL FUNCT (N^W^Fl) 

GO TO 170 
130 CONTINUE 
IEXIT-3 
GO TO 150 
l<tO CONTINUE 

IF ( I0IFF.EQ.2) GO TO 150 
IDIFF-2 
GO TO 50 
150 CONTINUE 

DO 160 I-1>N 
160 G( I)-M( IG-t-I) 

IF (IPRINT.EQ.O) RETURN 
PRINT 210# ITN# IFN# lEXIT 
PRINT 220# F 
PRINT 220# (X(I)#I>1#N) 

PRINT 220# (G(I)#I«1#N) 

RETURN 

170 CONTINUE 

IF (ITN.GE.l) RETURN 
CALL FUNCT (N#X#F) 

IFN-IFN+1 
CALL TWRIT 
DO 180 I=1#N 
Z-HH + XM( I ) 

ZZ»X( 1) 

X(I)«ZZ+Z 

CALL FUNCT (N#X#F1) 

W(IG+I)-(F1-F)/Z 
180 X( I ) -ZZ 
IFN*IFN+N 
GO TO 60 
C 

190 FORMAT (15H1ENTRY TO VAIOA#/) 

200 FORMAT ( 6H DFN »#F10.5#5H HH ■#F10.5#/#8H XM ( I ) «#(F10.5)) 
210 FORMAT (2AI5) 

220 FORMAT ( (8E15.7) ) 

230 FORMAT (1X#7HALPHA «#F10.5) 

END 
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SUBROUTINE MC IIB (A,N>IR) 

C OPTIMIZATION SUBROUTINE 
C FACTORIZE A MATRIX GIVEN IN A 
DIMENSION A(l) 

IR-N 

IF (N.GT.l) GO TO 

IF (A(l).GT.O.) RETURN 

A(1)«0. 

IR»0 
RETURN 
10 CONTINUE 
NP-N+1 
II-l 

DO 50 I»2/N 
AA-A( 1 1 ) 

NI»I I+NP-1 

IF (AA.GT.O.) GO TO 20 
A(II)=0. 

IR-IR-1 
II-NI+1 
GO TO 50 
20 CONTINUE 
IP-II+1 
II-NI+1 
JK-II 

DO ^0 IJ-IP/NI 
V»A( IJ ) / AA 
DO 30 IK-IJ,NI 
A( JK)-A( JK)-A(IK)*V 
30 JK-JK+1 
40 A(IJ )-V 
50 CONTINUE 

IF (A( II ) .GT.O. ) RETURN 
A( II ) -0. 

IR-IR-1 

RETURN 

END 


SUBROUTINE MCllE (A,N,Z#W/IR) 

C OPTIMIZATION SUBROUTINE 

C MULTIPLY A VECTOR Z BY THE INVERSE OF THE FACTORS GIVEN IN A 
DIMENSION A(l)> Z(l), W(l) 

IF (IR.LT.N) RETURN 
W ( 1 ) « Z ( 1 ) 

IF (N.GT.l) GO TO 10 
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Z(1)«Z(1)/A(1) 
RETURN 
10 CONTINUE 
DO 30 I»2>N 
U«I 
Il-I-l 
V»Z(I) 

DO 20 

V»V-A( IJ )*Z(J ) 
20 IJ-IJ+N-J 
W( I> -V 
30 Z(I)-V 

Z(N) *Z(N)/A(IJ ) 
NP-N+1 

DO 50 NIP»2>N 

I- NP-NIP 

II- IJ-NIP 
V«Z( 1)/A( II) 
IP-I+1 

IJ-I I 

DO <tO J»IP^N 
II-II+l 

40 V-V-A( II)*Z(J ) 
50 Z(I)«V 
RETURN 
END 


SUBROUTINE REFIN (ALFO) 

C HALVES MESH SIZE 

COMMON G(129,12>15)^S0(129>15)/E0(13i)^Z0(131),IV(129,15)#ITEl(15) 
ITE2(15)>A0(129)#Al(i29),A2(129),A3(129),B0(12)#Bl(12)j B2(12),B3( 
212)#Z(15),Cl(15),C2{15)^C3a5),XC(15)>XZ(15)#XZZ(15),YC(15),YZ(15) 
3#YZZ(15)#NX>NY#NZ#KT£1^KTE2>ISYM>K5YM,SCAL#SCALZ/YAW#CYAW^SYAW>ALP 
4HA>CAf S A, F MACH# Nl^N2iN3#I OUNCES #TSTEP> EPS 1>Q PR E( 129# 15 ) # SOPR E ( 129, 
515)#NQSTA#ZQSTA{15)#PCU1( 15)#PCU2(15)#PCQ3(15)#QQ1{15)#UQ2{15)#QQ3 
6(15)#QQ4(15)#PCS1(15)#PCS2(15 )#DSURF(15)#RDQ#RDS0#F00#FC1#F10#F11# 
7NDQ# IQ#KQ# AW ING# VOL DRG# I ORGP L T ( 129# 1 5 ) # SEC DRG ( 1 5 ) 

DIMENSION ALFO(l) 

MX-NX+1 

KY-NY+1 

MY-NY+2 

MZ-NZ+3 

MXO-NX/2+1 

MZO-NZ/2+2 

K»1 


79 



KK«K 

IF (KSYM.EQ.O) GO TO 10 

MZ0»NZ/2+3 

K«2 

KK»K 

10 DO 70 K-KK/MZO 
J-NY/2+1 
J J-KY 
20 I-MXO 
II-MX 

30 G(II>JJiK)»G( I#J#K) 

I- I-l 

II- I1-2 

IF ( I .GT.O) GO TO 30 
J*J-1 
J J-J J-2 

IF (J.GT.O) 60 TO 20 
DO 40 J«1/KY,2 
DO 40 I«2>NX#2 

40 G(I#J>K)».5*(G(I + 1# J#K)+G(I-1>J#K) ) 

DO 60 I-1#MX 
DO 50 J»2#NY>2 

50 G{I,J,K)«.5*(G(I,J+1/K)+G(I,J-1,K)) 

60 G(I,MY>K)»0. 

70 CONTINUE 

IF (KSYM.NE.O) GO TO 80 
MZM-MZO 
MZST«NZ+1 
GO TO 90 
80 MZM»MZ0 
MZST-MZ 
90 CONTINUE 

DO 100 J»1,MY 
DO 100 I*1»MX 

100 G( I, J#MZST) «G ( li J>MZM) 

IF (MZST.EO.l) GO TO 120 
MZST-MZST-1 
DO 110 J=1#MY 
DO 110 I-1>MX 

110 G( I,J#MZST)«0.5*(G( I^ J^MZM)+G(I,J,MZM-1 ) ) 
MZM=hZM-l 
MZST-MZST-l 
GO TO 90 
120 CONTINUE 

TYAW = SYAW/C YAW 
S1-. 5*SCAL 
N0«KTE1-1 
EO(NO) =0. 

K»2 

KK-K 

IF (KSYM.NE.O) KK=K+1 
DO 200 K«KK,MZ 
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N«NO 

I-MXO+1 

IF (K.LT.KTE1.0R.K.GT.KTE2) GO TO 150 
I1«ITE1(K) 

I2-ITE2(K) 

DO 130 1>11^12 
DSI-S0(I+1^K)-S0(I-1jK) 

DSK-SOt I>K+1)-S0( I>K-l> 

SX-AK I)*DSI 
SZ-CKK )*DSK 

DGI»G( I+1>KY#K)-G(I-1^KY>K) 

DGK«G( I»KY^K+1)-G(I»KY^K-1) 

R-AMINO( 1« IV( I«K) ) 

A-1.0-R+A0( I)*AO( I)+SO( I#K)*SO( I^K) 

H»R/ A 
FH«R*A 

AZ»-A0( I)*XZ(K)-SO(I>K)*YZ(K) 

BZ— A0( I)+YZ(K)+SO(I#K)*XZ(K) 

HZ-AZ*SX“BZ+FH*SZ 
FYY-l.O+SX+SX+H+HZ+HZ 
FXY-SX+H^ AZ+HZ 

U«A1(I)*DGI+CA*A0{I )+Sa*S0(I>K) 

W«C1 (K)+DGK+S YAW+CA*XZ(K)+SA+YZ(K) 

V«SA*A0( I)-CA*S0( I,K) 

130 6(I#KY+1jK)* 6{ I>KY-1#K)+(V*(1.0-H»BZ»HZ)-0*FXY-W*HZ )/ (FYY^flKKY) ) 
NO-NO+1 

EO(NO)«G( I2#KY#K)-G( Il#KYiK) 

N»NO 

I-Il 

IF (K.NE.KTE2.0R.YAW.LE.0. ) GO TO 150 
lAO I-I+l 

M-NX+2-I 

NO-NO+1 

EO(NO)«G(M>KY#K)-G(I^KY#K) 

IF (I.LT.MXO) GO TO lAO 
1*11 

150 I-I-l 

E-0. 

IF (IV( I,K) .NE.l) GO TO IBO 
Z2»Z(K)-TYAW*( XC (K) +S1*A0( I )*A0( I ) ) 

160 IF ( ZZ.GE.Z0< N-1) ) GO TO 170 
N«N-1 
GO TO 160 

170 R«(ZZ-Z0(N-1) )/(ZO(N)-ZO(N-l) ) 

E*R+EO(N )+{ l.-R )*E0(N-1) 

180 M-NX+2-I 

G(I>KY+1>K)«G(M>KY-1>K)-E 
G(M^KY+1,K) «G( I/KY-1,K)+E 
IF t IV< I,K) .NE.-l) GO TO 190 

G( I#KY#K)».5*G( li KY#K-l)+.25*{6( I#KY»K+1)+G(M>KY/K+1) ) 

IF (IV(I#K+l).LT.l) G( I#KY^K)«.5 *G{IjKY,K+ 1)+.25+(G(I,KY#K-1)+G(M> 
IKY^K-l) ) 
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G<M>KY>K) «G( I#KY>K) 

G( I»KY“1#K)». 5+(G(I,KY/K)+G(I#K Y-2» K ) ) 
G(«^ KY-1>K J «. 5+(G(M^ KY>K)+6(M>KY-2>K) ) 
190 IF ( I.GT.2) GO TO 150 
200 CONTINUE 

E0(N0+1)«0. 

K2-KTE1+(KTE2-KT£1) /2 
ALF0{KTE2)»ALF0(K2) 

K»K2-1 

KK-KTE2-1 

210 ALF0(KK)«.5*( ALFO(K )+ALFO(K+l) ) 
ALFO(KK-l) »ALFO(K ) 

KK-KK-2 

K-K-1 

IF (K.GE.KTEl) GO TO 210 

RETURN 

END 


SUBROUTINE SPLIF ( M, N^ S j F , FP, FP P , F P PP > KM> VM » KN> V Nf MOD E > F QM > IND ) 
C CUBIC SPLINE FIT WITH PRESCRIBED END CONDITIONS 
C INTEGRAL PLACED IN FPPP IF MODE GREATER THAN 0 

C IND SET TO ZERO IF DATA ILLEGAL 

DIMENSION S(l)i F(l), FP(1)» FPP(l), FPPP(l) 

I ND«0 

K«IABS(N-M) 

IF (K-1) 180^160#10 
10 K«(N-M)/K 
I-M 
J-M + K 

DS»S(J)-S(I) 

D-DS 

IF (DS) 20^180/20 
20 DF-(F( J )-F ( I) )/DS 
IF (KM-2) 30/40^50 
30 U».5 

V-3.*(0F-VM)/DS 
GO TO 80 
40 U»0. 

V-VM 

GO TO 80 
50 U— 1. 

V=-DS*VM 
GO TO 80 
60 I»J 
J-J+K 
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DS«S( J )-S( I) 

IF (D+DS) 180#180>70 
70 DF«(F(J )-F(I) )/DS 
B«1./(DS+0S+U) 

U»B*0S 

V«B+(6.+0F-V) 

80 FP(I)»U 
FPP{ I)«V 
U»(2.-U)+DS 
V«6.*DF+DS*V 
IF (J-N) 60>90#60 
90 IF (KN-2) 100#110,120 
100 V»(6.*VN-V) /U 
GO TO 130 
110 V-VN 

GO TO 130 

120 V»(DS + VN + FPP{ I) )/(l.+FP(I) ) 

130 B«V 
O-OS 

140 DS»S(J)-S( I) 

U*FPP(I)-FP(I)*V 
FPPP( I) »(V-U ) /DS 
FPP( I)«U 

FP(I ) *( F ( J )-F { I ) ) /DS-DS*{ V+U+U) /6. 

V«U 
J»I 
I =I-K 

IF (J-M) 140>150,140 
150 I-N-K 

FPPP(N)»FPPP( I) 

FPP(N)=B 

FP{N)«DF+D*(FPP(I)+6+6)/6. 

IND-1 

IF (MODE) ia0»180#160 
160 FPPP(J)-FQM 
V-FPP ( J ) 

170 I«J 
J »J + K 

DS«S( J)-S( I) 

U-FPP ( J ) 

FPPP( J) »FPPP{ I ) + .5*DS+(F( I )+F( J)“03*0S*(U + V ) /12. ) 
V-U 

IF (J-N) 170»180#170 
180 RETURN 
END 



SUBROUTINE INTPL ( M I # N I > S I> F I > M/ N# F # F P, F R P > FP P MODE ) 

C INTERPOLATION OF CUBIC SPLINE BY TAYLOR SERIES 
C ADDS CORRECTION FOR PIECEWISE CONSTANT FOURTH DERIVIATIVE 
C IF MODE GREATER THAN 0 

DIMENSION SKI)/ FI(1)> S(l)/ F(l)/ FP(1)/ FPP(l)/ FPPP(l) 
K«IABS(N-M) 

K-(N-M)/K 

I«M 

MIN-MI 
NIN-NI 
D«S (N)“S (M) 

IF (D+(SI(NI)-SI(MI))) 10/20/20 
10 MIN-NI 
N1N»MI 

20 KI-IABS (NIN-MIN) 

IF (KI) 40/40/30 
30 KI»(NIN-MIN)/KI 
40 II-MIN-KI 
C*0. 

IF (MODE) 60/60/50 
50 C«l. 

60 Il-II+KI 
SS-SI (II ) 

70 I-I+K 

IF (1-N) 80/90/80 
80 IF (D*(S(I)-SS) ) 70/70/90 
90 J»I 
I-I-K 

SS«SS-S ( I ) 

FPPPP»C*(FPPP(J)-FPPP(I))/(S(J)“S(D) 

FF«FPPP(I)+.25*SS*FPPPP 

FF-FPP ( I)+SS*FF/3. 

FF»FP( I )+.5*SS*FF 
FI(II)»F(I)+SS*FF 
IF (II-NIN) 60/100/60 
100 RETURN 
END 


SUBROUTINE THREED ( I P LOT / S V/ SM/ C P/ X / Y/ T I TL E / Y A/ AL/ V LD/ C L/ C D/ CHDRDO 
1/XSCAL/PSCAL/ LABEL/NIT/UC/VC/WC/NFl) 

C GENERATES THREE DIMENSIONAL PLOTS 
C GENERATES CALCOMP PLOTS ON CDC 6600 

COMMON G(129/12/15)/S0(129/15)/E0(131)/Z0( 131) / IV(129/15) / 1TEK15) 
1/ ITE2( 15)/ A0( 129)/ AK129)/ A2(129)/A3(129)/60(12 )/B1 (12)/B2(12)/B3( 
212)/Z(15)/C1(15)/C2(15)/C3(15)/XC(15)/XZ(15)/XZZ(15)/YC(15)/Y2(15) 
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3^YZZ(15)>NX#NY>NZ>KTEl,KT£2>lSYM#KSYM,SCAL,SCALZ#YAW#CYAW#i)YAw^ALP 
4HA,CA#SA>FMACH>N1#N2#N3# 10#N0ES/TSTEP>EPS1>QPR£( 129>15)# SOPRt ( 129, 
515),NQSTA,ZQSTA( 15),PCQi(15),PCQ2(15),PCQ3(15),QQl( 15 ),QQ2 (15),UQ3 
6(15),QQ4(15)#PCSl(15),PCS2{15),0SURF(15),R0Q>RDS0,F0O,F01, F10,F11, 
7N0Q, IQ,KQ, AWING, VQLDRG, I ORGP L T ( 129, 15 ) , SEC DRG ( 1 5 ) 

DIMENSION X(l), Y(l), SV(1), SM(1), CP(1), TITLE(IO), R(20), UC(1) 
1, VC(1), WC(1) 

M»1 

LX-NX/2+1 

MX«NX+1 

MY-NY+2 

IF ( XSCAL.NE.O. ) SCALX- . 5* ABS ( XSCAL ) /CHQRDO 
IF (PSCAL .GE. 0. ) SCALX»5./(Z(KTE2)-Z(KTE1) ) 

SCALP — 1.00 

IF (PSCAL.NE.O. ) SCALP--. 5/ABS( PSCAL ) 

TX-3.0 

SX — SCALX*XC(KTE1) 

IF (IPLQT.NE.l) GO TO 10 

CALL PL0TS8L ( 10000, 22HJ E F F MCFADDEN 3210WWH) 

10 IPLOT-0 
CALL FRAME 

CALL PLOT (1.25,1., -3) 

ASRAT»2.+(Z(KT£2)-Z (KTEl) )**2/AWING 
ENCODE (60,90,R) FM ACH, C L , VOL DR G, ASR AT 
CALL SYMBOL ( . 50, 0 . 75, . 14, R , 0 . , 60 ) 

ENCODE (60, 100, R) 

CALL SYMBOL ( . 50, 1 . 2 5, . 14, R, 0 . , 60 ) 

20 CONTINUE 
K-1 

IF (KTE2.LT,10) K-2 
30 K-K+2 

IF (KTE2.lt. 10) K«K-1 
IF (K.GT.KTE2) GO TO 70 
IF (K.LT.KTEl) GO TO 30 

11- ITEKK) 

12- ITE2(K) 

CALL VELO (K,K,SV,SM,CP,X, Y,UC, VC,wC ) 
SY-5.*(Z(K)-Z(KTEl))/(Z(KTE2)-Z(KTEl))+2.45 
SCP«5.*(Z(K)-Z(KTEl))/(Z(KTE2)-Z(KTEl))+2.75 
DO 40 1-11,12 
X(I)-SCALX*X( I )+SX 
Y(I)-SCALX*Y(I)+SY 
40 CP(I)«SCALP*CP(I)+SCP 
IF (M.EQ.2) GO TO 50 
N-I2-LX+1 

CALL LINE (X(LX),CP(LX),N,1,0,1,0.,1.,0.,1.) 

GO TO 30 
50 N-I2-I1+1 

DO 60 1-11,12 
60 X(I)«X(I)+TX 
N-I2-I1+1 

CALL LINE (X( I1),Y( I1),N,1,0,1,0.,1.,0.,1. ) 
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GO TO 30 
70 CONTINOE 
M*M+1 

IF (M.GT.2) GO TO 80 
60 TO 20 

80 CALL DRAGC (1#SCALX) 

I0«1 

CALL PLOT (-1.25>-l.#-3) 

RETURN 

C 

90 FORMAT {-itHM ■ # F3 . 2# IH^ 2X ^ 5HCL » > F 3 . 2# 2X^ 6HC DW « # F5 . 4, IH# 2X> A 

IHA * »F3.1) 

100 FORMAT (22HUPPER SURFACE PRESSURE, 5X, 15HWING AND SHOCKS) 

END 


SUBROUTINE RE ADOS (NQSTA, ZQSTA,PCQ1, PC0 2, PCQ3, 0 1, 02, 0 3, QA, PC S 1, PC S 
12,DSURF,FMACH) 

C READS IN ASSIGNED SPEED DISTRIBUTION 

DIMENSION ZOSTA(l), PCOl(l), PC02(1), PC03(1), 01(1), 02(1), 03(1) 
1, 04(1), PCSKl), PCS2(1), DSURF(l) 

IREAD*9 

IWRIT*6 

AA0-1./FMACH*=*‘2+.2 
WRITE (IWPIT,20) 

READ (IREAD,4C) 

READ (IR£AD,40) 

READ (IREAD,50) FOSTA 
NOSTA-FQSTA 
DO 10 K*1,N0STA 
READ (IR£AD,40) 

READ (IREAD,50) Z0STA(K),PC01(K),PC02(K),PCQ3(K),01(K),02(K),Q3(K) 
1,04(K ) 

READ (IREAD,40) 

READ (IR£AD,50) PCS 1 ( K ) , PCS2 ( K ) , DSUR F ( K ) 

WRITE (I WRIT, 30) K,ZQSTA(K),PCQ1(K),PCQ2(K),PCQ3(K),Q1(K),Q2(K),Q3 
1(K),04(K),PCS1(K),PCS2(K),DSURF(K) 

Q1(K)-SQRT( (AA0*Qi(K)**2) /(!. + . 2*01 (K)* + 2) ) 

02(K)«S0RT( (AA0*Q2(K)**2)/(1.+.2*02(K)**2) ) 

Q3(K)«SQRT( (AA0*Q3(K)**2)/ ( 1 . + . 2*0 3 ( K ) ♦* 2 ) ) 

Q4(K)«SQRT( (AA0*04(K)**2) / ( 1 . +. 2*04 ( K )**2 ) ) 

10 CONTINUE 
RETURN 
C 

20 FORMAT (55H1 PARAMETERS TO DEFINE THE ASSIGNED DESIGN MACH NUMBER 
113HDISTRIBUTI0N., //,3H K, 9X, 1HZ,6X, 4HPCM1 , 6X, 4 HPCM2, 6X, 4HPCM3> 8X, 
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22HM1^ 8X>2HM2> 8X> 2HM3^ 8X> 2HM^, 6X> 4HPC XI# 6X# ^HPCX2# 5X# 5HDSURP# / ) 
30 FORMAT ( IX # 12 # IIF 10. 5 ) 

40 FORMAT (IX) 

50 FORMAT (8F10.5) 

END 


SUBROUTINE SETQS ( N E # N X# Q PRE # SO# SOP RE # ITEl# ITE2#KTE1#KTE2# Z#ZQSTA# 
1A0#PCQ1#PCQ2# PCQ3#SI#FI#Q1#Q2#Q3#Q4#PCS1#PCS2#DSURF#NQSTA) 

C DEFINES ASSIGNED SPEED DISTRIBUTION BY EXPONENTIAL SPLINE 
C ALLOWS EASY CONSTRUCTION OF SHOCKLESS DISTRIBUTION USING t AND G 
DIMENSION QPRE(NE#D# SO(NE#D# SOPREtNE#!)# ITEKD# ITE2U)# Z(1 
1)# ZQSTA(l)# AO(D# PCQKD# PCQ2(D# PCQ3(D# SKI)# FKD# Ql(l) 
2# U2(D# Q3(D# Q4(l)# PCSKD# PCS2(D# DSURF(l) 

DIMENSION X(4)# Y(4)# E(4)# G(4)# A(4)# B(4)# C(4)# D(4) 

DATA ISET/O/ 

BUMP(X)-16.^{X*(l.-X))++2 

E ( 1) « .007 

E(2)a.007 

E ( 3) -.007 

G(1)«0. 

G(2)«-3. 

G( 3) «45. 

LX-NX/2+1 

MX-NX+1 

KM»2 

VM-0. 

KN*3 

VN-0. 

K2-2 

K*KTE1-1 
10 K«K+1 

Il-ITEl (K) 

I2»ITE2(K) 

ZO-Z(K) 

K2-K2-1 
20 K2-K2+1 
K1-K2-1 
Z1-ZQSTA(K1 ) 

Z2»Z0STA(K2) 

IF (Z0.GT.Z2.AND.K2.LT.NQSTA) GO TO 20 
R1«(Z2-Z0)/(Z2-Z1) 

R2-1.-R1 

DLEN=AO( I2)-A0(LX) 

PXl-Rl+PCQl (K1)+R2*PC01(K2) 

PX2»R1*PCQ2(K1)+R2*PCQ2(K2) 
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PX3»R1«PCQ3(K1)+R2*PCQ3{K2) 
X(1)«A0(LX)+PX1*0LEN 
X(2)*A0(LX)+PX2+DLEN 
X(3)»A0( LX)+PX3*DLEN 
X(^)«A0(I2) 

Y(1)«R1*U1(K1)+R2*Q1(K2) 

Y(2)«R1*Q2 (K1 )+R2+Q2(K2) 

Y(3)«Ri+C3{Kl )+R2*Q3(K2) 

Y(<f)*Rl*Q<t(Kl)+R2 + Q^(K2) 

CALL SPTEN ( A # X# Y^ E ^ G> A^ B> C# 0,KM^ VM# KN, VN ) 
DO 30 1=1, MX 
QPRE( I,K)«0. 

IF (ISET.LQ.O) SOPRE( I,K)-SO(I,K) 

30 CONTINUE 
LS = I1 

AO LS-LS+1 

IF (AO(LS).LT.X(D) GO TO AO 
NN-MX-LS+1 
DO 50 1=1, NN 
J-LS+I-1 
50 Sl( I)«AO(J ) 

CALL INTEN ( NN , S I , F I , A , X , A , B , C, 0 , E, G ) 

DO 60 I=LSiMX 
J*I-LS+1 

60 QPRE( I,K)»FI( J ) 

DENO-1, /FLOAT (LS-Il ) 

DO 70 I«IX,LS 

70 QPRE { I,K) -FLOAT ( I-I 1 ) ♦ DEN0+ Y ( 1) 

IF (ISET.tQ.l) GO TO 90 

PX2»R1*PCS2{K1)+R2*PCS2(K2) 

PX1«R1*PCS1(K1)+R2*PCS1(K2) 

DS0PR£«R1*DSURF(K1 )+R2+DSURF(K2 ) 

X1=A0(LX)+PX1*DLEN 

X2=A0(LX)+PX2*DLEN 

I-Il 

80 1=1+1 

IF ( A0( I ) ,GT. X2) GO TO 90 
IF (AO(I).LT.Xl) GO TO 80 
XX = ( A0( D-Xl) /(X2-X1) 

SOPRE ( I,K)«SC{ I,K)-DSOPRE*BUMP( XX) 

GO TO 80 

90 IF (K.LT.KTE2) GO TO 10 
ISET-1 
RETURN 
END 
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SUBROUTINE SPTEN ( N^ S> F # E> 6# A> 6>C KM» VM> KN> VN ) 

C COMPUTES EXPONENTIAL SPLINE WITH WEIGHTING FACTORS 
C E IS TENSION PARAMETER# SMALLER E PRODUCES LESS OSCILLATION 

C G IS WEIGHT FACTOR# LARGER G PRODUCES MORE SAG 

DIMENSION S(D# F(D# Ed)# G(D# A(D# B(l)# Cd)# Dd) 

NM-N-1 

H-S(2)-Sd) 

HI-l./H 

SI-l./SINH(H/Ed) ) 

TI-SI»COSH(H/Ed) ) 

IF (KM-2) 10#20#30 
10 Bd)»Ed)+(E( D+HI-TI) 

Cd)-Ed) + (SI-Ed)*HI) 

Dd)»VM + HI*(F d)-F(2) )+E(l)*Gd )*(SI-TI) + .5 + H*Gd) 

GO TO 40 
20 Bd)>l. 

C d)*0. 

Dd) »VM 
GU TO 40 
30 B(1)»-TI 
C d) »SI 

Dd)«Gd)*(SI-TI)+VM*Ed) 

40 XX-l./Bd) 

Dd)«XX*D(l) 

D0 50 I»2#NM 

HM«H 

HIM*HI 

SIM«SI 

TIM«TI 

H*Sd + l)-Sd) 

HI-l./H 

SI-l./SINH(H/Ed) ) 

TI-SI*CQSH{H/E(I ) ) 

A( I) -E { I-l )♦( SIM-HIM+E ( I-l ) ) 

Bd)«E( I)*(HI^Ed )-7I)+£( I-1)*{HIM*E ( I-D-TIM) 

C ( I ) *E ( I )*(SI-HI*E( I ) ) 

Dd) *HIM»( Fd )-Fd-l ) )+E( I-1)*G( I-l ) ♦ ( S IM-T IM ) + . 5 + HM* G ( I-l ) +h I* ( F { 
II )-F( I + l ) )+E ( !)♦& (I )*(S I-TI ) + .5*H*G(I ) 

Cd-l)»XX*Cd-l) 

XX-l./(B(I)-A(I)*Cd-l) ) 

50 0(I)»(D(I)-A( I)*D(I-1))’*‘XX 
IF (KN-2) 60#70#80 
60 A(N)-E(NM)+(HI*E(NM)-SI ) 

B{N)»E(NM)*(T I-HI*E(NM) ) 

D(N)«VN + HI+(F (NM)-F(N) ) -. 5*H* G( NM ) +E { NM ) *G ( NM ) * ( T I-S I ) 

GO TO 90 
70 A(N)«0. 

B(N) »1. 

D(N) »VN 
60 TO 90 

80 A(N)»- SI/E(NM) 

B(N)=TI/E(NM) 


89 



D(N)=VN+G(NM)*(TI-SI )/E(NM) 

90 C(NM) »XX*C (NM ) 

XX»1./(B(N)-A(N) + C(NM) ) 

D(N>«(D(N)-A(N)*D(NM) )*XX 

DO 100 J«2>N 

I«N+1-J 

100 D(I)=D(I)-C(I)*D(I+1) 

DI»D(1) 

DO 110 I»1>NM 
DIP»D(I+1) 

H-S( I+1)“S( I) 

SI«1./SINH(H/E(I ) ) 

A(I )=F(I+1)-F (I)+(DI-DIP)*E<I )**2-.5*G(I)*H*+2 
BtI)*F(I)+(G(I)-DI)*E(I)+*2 
C{I)«SI*(DIP“G(I) )*E(I)**2 
D(I)=Sl*(DI-G(I))*E(I)+*2 

110 DI-DIP 
RETURN 
END 


SUBROUTINE INTEN (NX# SI^ F I#N, S# A, B,C^ D> E> G) 

C COMPUTES INTERPOLATED VALUES FROM EXPONENTIAL SPLINE 

DIMENSION SKI)# FKD# S(D# A(D# B(D# C(l)# D(D# E(D# 6(1) 

I-O 

J-1 

JP-J+1 
10 I-I+l 

20 IF (SI(I).LT.S(JP)) GO TO 30 
IF (JP.EQ.N) GO TO 30 
J -J + 1 
JP-J+1 
GO TO 20 

30 HI-1./(S(JP)-S(J) ) 

EI-1./E(J) 

FI{I)-B(J)+HI*A(J)*(SI(I)-S(J))+C(J)*SINH{E1*(SI(I)-S(J)))+D(J)*SI 
1NH(EI*(S( JP)-SI(I)) )+.5*G(J)«(SI(I)-S(J) )+*2 
IF (I.LT.NX) GO TO 10 
RETURN 
END 
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C THIS VERSION OF SUBROUTINE YSWEEP IS VECTORIZED AND CAN BE INVOKED 
C BY SETTING FSWEEP TO -1.0 ON THE APPROPRIATE INPUT CARD 
SUBROUTINE VYSWEEP 
C ROW RELAXATION 

COMMON G(129# 12> 15 ) # S0( 129, 1& ) , E0( 131 ) # Z0( 1 31 ) # I V ( 129> 1 5 ) / IT 1 1 ( 15 ) 
1, ITE2(15),A0(129),A1(129), A2(129),A3(129),B0{12),B1(12), B2(12),B3( 
212),Z(15),C1(15),C2(15),C3(15),XC(15),XZ{15),XZZ(15),YC(15),YZ(15) 
3, YZZ( 15) ,NX,NY,NZ,KTE1,KTE2,ISYM,KSYM,SCAL,SCALZ,YAW,CYAW, S YAW, ALP 
^HA,CA,SA,FMACH,Ni,N2,N3, I0,ND£S,TSTEP,EPS1,QPR£(129,15),S0PRE(129, 
515),NQSTA,2QSTA(15),PCQ1<15),PCQ2(15),PCQ3( 1 5 ) , QQl ( 15 ) , QU2 ( 15 ) , Q03 
6{15),QQ4(15),PCSi(15),PCS2<15),DSURF(15),RDQ,RDS0,F00,F01, F10,F11, 
7NDQ,IQ,KQ,AWING, VOLDRG,IDRGPLT( 129,15),SECDRG(15) 

COMMON /FLO/ STR I P, P 1, P 2, P 3, B E T A, FR , 1 R, J R, K R, DG, I G, J G, KG, N S 
COMMON /SWP/ DXYZ(129),GK1(129,15),GK2(129,15),SX(129),SZ(1Z9),SXX 
1(129),SXZ(129),SZZ(129),R0(129),R1(129),C(129),D(129),G10(15),620( 
215),G30(15),G40(15) ,G1( 15),G2(15),I1, I2,K,L,N0, LX,MX,KY,MY,T1, AAO, 
3Q1,Q2,TYAW,S1 

COMMON /DIM/ NX1,NY1,NZ1, FDIM 
COMMON /CRAY/ KA 

COMMON /VECT/ YP(129),3AVE(129),TEMP2(129),TEMP(129),TEMP1(129),AJ 
1(129),H(129),FH(129),aZ( 129),BZ(129),CZ(129),DZ(129),DGI{129),DGJ( 
2129),DGK(129),U(129),V(129),W(129),AU(129),AV(129),Q0(129),HGLD1{1 
329),H0L02{129),H0LD3(129),HQLD4(129),H0LD5(129),hZ(129),AA(129),FX 
9X(129),FYY(129),FXY(129),BV(129),UU(129),VV(129),WW(129),UV(129),U 
5W(129),VW(129),AZZ(129),AXX(129),AXZ<129),R(129),AXTU29),AYT(129) 
6,AZT(129),DGII(129),DGJJ(129),DGIJ(129),DGIK(129),DGJK(129),AC(129 

7) ,AB(129),AYY(129),AYZ( 129 ) , BP ( 129 ) , BM ( 129),B(129),AXY(129),CG(129 

8) ,DGKK(129),A(129),S(129) 

BETX-,01 

BETY-.15 
BETZ*. 1 

BSCAL*!./ (l.+FDIM) 

BSCALl*!./ (2.*(1. + FDIM) ) 

Jl«2 

IF ( FMACH.GE.l. ) Jl«3 
C ( Il-l) «0. 

D( Il-l) *0. 

DO 10 1-11,12 
R0( I)-l. 

RKI ) -1. 

GK1(I,1)-G(I,1,L) 

10 GKK I,J1-1)-G(I,J1-1,L) 

J-Jl 

13-12 

20 BC — Tl+BK J )*C1(K ) 

DO 30 1-11,13 
YP(I)«SO(I,K)+BO(J) 


91 


SAVE ( I ) -1 .0-R0( I ) 

TEMP2(I)«YP(I)*YP(I) 

TEMP(I)»AC(I)*AO(I) 

AJ(I)«SAVE( I)+TEMP2(I)+TEMP(I) 

30 CONTINUE 

DO 40 I»I1#I3 
H(I)-R0(I)/AJ(1) 

FH(I)-RO(I)+AJ(I) 

TEMP1(I)«A0(I )*(4.*TEMP2(I)-FH( I) ) 
TEMP2(I)*YP(I)*(4.+TEMP(I)-FH(I)) 

AT«XZ(K)+XZ(K )-YZ(K)+YZ(K ) 

BT»(XZ(K)+XZ(K))*YZ(K) 

AZ(I )*-A0( I )*XZ(K )-YP(I )*YZ<K) 

BZ( I ) — A0( I )+YZ (K)+YP ( I )+XZ (K) 

TEMP(I)»H(I)*H(I) 

CZ(I )=TtMP(I)*(TErPl(l )*AT-TEMP2(I)*BT)-A0( I)*XZZ(K)-YP(I)*YZZ(K) 
DZ( I)«TEMP (I)+(TEMP2( I)*AT+TEMP1{I)*BT)-A0( I)*YZZ(K)+YP(I)+XZZ(K) 
40 CONTINUE 

DO 50 I«I1>I3 

OGI(I)=G( I+l# J>U-G(I-1#J#L) 

DGJ( I)»G{I>J+1#L)-6K1(I#J-1) 

UGK( I)»G(I# Ji L + D-GKK I>J ) 

50 CONTINUE 

00 60 I=I1#I3 
TEMP1(I)«A1(I )*DGI(1) 

TEMP2(I)«-B1(J)*DGJ(I) 

U(I)*TtMPl{I)-SX(I) + TEMP2(I)+CA>»‘A0(I)+SA4YP(I) 

V(I)»TEMP2( I)+SA*AO( l)-CA^YP(I) 

W(I)«RO( I)'*‘(Ci(K)*DGK(I)-SZ(I)*T£MP2(I) + SYAW + CA*XZ(K)+SA4YZ(K)+H( I 
1)*(U(I)»AZ(I)+V(I)*BZ(I))) 

AU(I)“U(I)+W(I)«AZ(I) 

AV(I)«V(I)+W( I) + BZ(I) 

TEMP< I)«H(I)>MU(1)*U(I)+V(I)*V(I)) 

QQ(I )*TEMP(1)+W(I )*W(I ) 

60 CONTINUE 

DO 70 I»I1jI3 
HOLDK I)-.2*QQ( I) 

AA(I)»DIM(AA0>H0LD1(1) ) 

HZ(I)»AZ( I)«SX(I)-BZ(I)+FH(I)*SZ<I) 

FXX( I)-1.0 + H( I)*AZ( I)*AZ( I) 

FYY(I )»1.0+SX(I)*SX(I )+H(I )+HZ(I)*HZ(I) 

FXY( I)-SX(I)+H(I)*AZ(I)+HZ(I) 
BV(I)-AV(I)-AU(I)*SX(I)-FH(I)+Vi(I)*SZ(I) 

70 CONTINUE 

DO BC I=I1>I3 
UU(I)«H(I)tAU(I)4AU(I) 

VV(I)*H(I)*BV(I)+BV(I) 

WW(I )«FH(I)*W(I)=FW(I ) 

UV(I)«H{I)*AU(I)+BV(I) 

UW(I)>AU(I)*W(I) 

VW(I)»BV(I)*W(I) 

AXX( I)»R1(I)+(FXX(I)*AA(I)-UU(I)) 
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A2Z(I)«FH(I)+AA(I)-WW(I) 

AXZ( I)» (2.0*R0( I )♦( AZ( I )*AA( I )-UW(I ) ) ) 

80 CONTINUE 
DO 90 

H0LD1(I)»-TEMP2(I)*(AXX(I)*SXX(I)+AZZ(I)^SZZ(I)+AXZ(I)*SXZ(I) ) 
H0LD2( I) »( AA( !)♦ (CZ (l)*TEMPi(I)+(DZ(I)-SX(I)*CZ(I))*TEMP2(I)))*R0( 
II ) 

HQLD3( I ) »CA*( AU( I )*AU( I )-AV ( I ) ♦ AV ( I ) ) + ( S A+ S A ) ♦ AU ( I ) ♦ A V ( I ) 
HOLDA(I)»TEMP(I)+(U(I)*AO(I)+V(I)+YP(I)+2.0+W(I )♦( AO(I)+AZ(I)+YP( I 
1)*BZ( I) ) ) 

H0LD5(I)»-WW( I)+(CA+XZZ(K)+SA+YZZ<K))-W(I)*W(I)*(U(I)+CZ(I)+V(I)*0 
1Z(I) ) 

R( I) -HOLDl ( I) +T1*(H0LD2( I)-H( I )*(H0LD3( I)-HQLDA( I ) )+H0LD5( I) ) 

90 CONTINUE 

DO 100 1*11,13 
AXT(I )*AU{I)*A1(I ) 

AXT( I)»ABS(AXT(I) ) 

AYT( I)»BV(I)*B1(J ) 

AYT(I)*ABS(AYT(1) ) 

AZT( I)-FH( I )+W( I )*C1 (K) 

AZT( I)*ABS(AZT(I) ) 

SAVE(I)=AMAX1(AXT(I),AYT(I),AZT(I),{1.-R0(I))) 

HOLDK I ) -R0( I )*BETA+AA( I) / $AVE( I ) 

AXT(I)*AXT(I)*H0LD1(I) 

AYT(I)*AYT(I)*H0LD1(I) 

AZT(I)-AZT(I)*H0LD1(I) 

100 CONTINUE 

DO 110 1*11,13 

DGII(I)*G(I+1,J,L)-G(I,J,L)-G(I,J,L)+G(I-1,J,L)+A3(I)+DGI{I) 

DGJJ(I)-G(I,J+1,L)-G(I,J,L)-6<I,J,L>+G(I,J-1,L)-B3(J)*DGJ(I) 

DGKK( I)*G( I,J,L+1)-G(I, J,L)-G(I,J,L)+G<1,J,L-1)+C3(K)*DGK(I) 
DGIJ(I)=G(I+1,J+1,L)-G(I-1,J+1,L)-G(I+1,J-1,L)+G(I-1,J-1,L) 

DGIK( I ) -G( I + l, J, L+1 )-G( I + l, J, L-l)-G( 1 -1 , J, L + 1 ) +G ( I-l , J , L-1 ) 
0GJK(I)*G(I,J+1,L+1)-G(I,J-1,L+1)-G{I,J+1,L-1)+G{I,J-1,L-1) 

ACd ) *T1*A1(I )*C1(K) 

AB( I)--T1*A1( I)fBl(J) 

AXX(I)«AXX(I)^A2(I) 

AYYd )* (FYY( I )*AA (I )-VV(I))*B2(J) 

AZZ( I)*AZZd)+C2(K) 

AXY( I)*“R1(I)*(FXY(I) + AA( I)+UV( I))*(ABd) + ABd)) 
AXZd)»AXZd)*ACd) 

AYZ( I)*-RO( I)*(HZ(I)*AA( I)+VWd ) )*(BC+BC) 

BP ( I ) *AXX ( I ) 

BMd)*AXXd) 

B( I) *-AXX( I)-AXX(I)-Ql+(AYYd) + AZZd)) 

SAVE ( I)*AXXd )*DGII ( I)+AYYd)*DGJJ (I )+AZZ( I)»DGKK( I ) + AXY(I )*DGIJ ( I 
l)+AYZd )*DGJK ( I ) +AXZ ( I )*DGIKd ) 

110 CONTINUE 

DO 120 1-11,13 

IF (OQ(I).LT.AA(I)) GO TO 120 
NS-NS+1 

Sd)*SIGN(l.,Ud) ) 
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IM»I-IFIX(S{I ) ) 

IMM«IM^IFIX(S(D) 

AXX(I)«UU(I)+A2(I) 

AYY(I )=VV(I)+B2{J ) 

AZZ( I)«WW(I)+C2(K) 

AXY(1)«8.*S(I)*UV(I)*AB(I) 

AXZ(I)=8.*S(I )+UW(I )*AC(I ) 

AYZ( I)«8.*VW( I)*BC 

H0LD1(I)-(FXX{I)*QQ(I)-UU(I))*A2(I) 

H0LD2(I ) = (FYYd )*QQ (I )-VV(I ) )»B2( J ) 

H0LD3( I ) »(FH( I)>M3Q( I )-WW ( I ) )*C2 ( K ) 

HOLDA( I ) »-(FXY( I )*QQ( I) +UV( I) )*(AB{I)+AB(D) 

HDLD5(I)«(AZ( I)»QQ(I)-UW(I) )*{AC(I)+AC(I) ) 

TEMP ( I ) »-( HZ( I)*QQ( I )+VW( I ) )♦ (BC+BC ) 

TEMPK I)-AA(I )/QQ{I) 

T£MP2(I )»HaLDl(I )*DGII (I )+H0LD2( I)*DGJJ ( I ) +H0LD3 ( I ) *DGKK ( I ) +H0LD4 ( 
II )*DGIJ ( I )+TEMP( I )*DGJK( 1)+H0LD5 ( I )*DGIK( I ) 
DGII(I)«G(I#J^L)-G(IM»J#L)-G<IM#J#L)+6(IMM,J>L)+A3(I)ADGI(1) 
DGJJ(I)»G{I>J>L)-G( J-1# L)-G(I> J-1#L ) + GKl ( I, J-2 ) -B 3 ( J ) ♦ DG J (I) 
DGKK(I) =G( dJ#L)-G(I,J#L-l)-G(I# J,L“i)+GK2 ( I > J ) +C 3 ( K ) ♦DGK (I ) 

DGIJ{ I)-G( dJ>L)“G( IM#J/L)-G(I#J-1>L)+6(IM# J-1>L) 

OGIK( 1)«G( d J»L)-G( I>J>L-l)-G(IMf J#L)+G(IM> J>L-1) 

DGJKd )»G( d J»L)-G( I>J#L-1)“G(I, J-1 jL )+G(I> J- 1#L-1) 
TEMP(I)»AXX(I)*DGII{I)+AYY(I)+DGJJd)+AZZd)*DGKK(I) + AXYd)*DGlJ(I 
1 ) + AYZ ( I )*DGJK d ) + AXZd )*DG1K{ I) 
B(I)».5*{TEMPl{I)-l.)*{AXXd)+AXXd)+AXZd)+AXY(I)) 
BPd)*TEMPld)*HOLDl(I)-(l.-S(I))*B(I) 

BM( I) -TEMPld )*HQLD1( I)-{ 1.0+S( I) )*B( I ) 

B( I) --TEMP Id )*(H0LD1( I )+HULDl( I) +Q2* ( H0LD2 ( I ) +HOLD 3 { I ) ) ) + (TEMPl( I 

I) -1. )♦{ 2. + ( AXXd )+AYY(I ) + AZZ( I) )+AXY( I ) +AYZ ( I ) + AXZ( I ) ) 

SAVE( I )»(TEMP1( 1 )-l. )*TEMP( I)+TEMP1( I )*TEKP2( I) 

120 CONTINUE 

DO 130 I-11#I3 
130 Rd ) -R ( I ) + SAVE (I ) 

DO lAO I»Ild3 

IF (ABS(R(D) .LE.ABS(FR)) GO TO lAO 
FR-Rd) 

IR-I 
JR » J 
KR-K 

140 CONTINUE 

DO 150 I-Il>13 

Rd)»Rd)-AYTd)*(GKld^J-l)-Gd>J-l#L))-AZTd)=«c(GKld>J)-Gd>J>L- 

II) ) 

Bd)»Bd )-AXT d)-AYT( I)-AZT(I) 

BMd)»BMd)+AXT(l) 

150 CONTINUE 

DO 160 1-11,13 

Bd) «1.0/(B(I )-BMd )*Cd-l) ) 

Cd) -Bd ) + BPd) 

160 Od) -8(1 ) + (Rd)-BMd )*Dd-l) ) 

I-I3 
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CG(I3+1)»0. 

DO 170 M-Il>13 

CG(I )«D( I )“C( I)*CG( I+l) 

GK2( I>J )>GK1( I>J) 

GKK I#J )«G( J>L) 

G(l^J^L)«G(Ii JfD'CGd) 

170 I-I-l 
I-I3 

DQ 180 M«I1>I3 

IF ( ABS(CG(I) ).LE.A8S(DG) ) GO TO 180 
DG»CG( I ) 

IG»I 

JG-J 

KG«K 

180 I-I-l 
J«J + 1 

IF (J-KY) 20/190#210 
190 IF ( 12.GT. ITE2(K) ) I3«ITE2(K) 

IF ( ITE2(K ) .EQ.MX ) I3-LX 
DO 200 I«I1>I3 
LV«IABS(1-IABS( IV(I>K) ) ) 

RO(I)»AMINO(LV>IABS(IV(I>K))) 

200 RKI )*LV 
GO TO 20 
210 N-NO 
I-LX+1 

IF (K.LT.KTE1.0R.K.GT.KTE2) GO TO 230 

I0-NX+2-I3 

DO 220 I»I0fI3 

AJ ( I )«1.“R0( I )+A0{ I ) + A0( I )+S0 (I/K)*SO(I,K) 

H( I ) =R0( I ) /AJ (I) 

FH(I)-RO(I)*AJ(I) 

AZ(I)»-AO(I)+XZ(K)-SO(I#K)*YZ(K) 

BZ(I )«-A0( I)*YZ(K)+SO(I#K)*XZ(K) 
HZ(I)»AZ(I)+SX(I)-BZ(I)+FH(1)+SZ(I) 

FYY(I)«1. + SX( I)*SX( I )+H(I )+HZ(I )^HZ< I ) 
FXY(I)»SX(I)+H{I)+AZ(I)*HZ(I) 

DGK I ) -G( I + 1,KY>L )-G( I-1,KYjL ) 

DGK( I)«G( I#KY#L+1)-GK2( I#KY) 

V(I)»SA*AO(I)-CA*SO(IfK) 

U(I)-A1( I)*DGI(I)+CA*AO( I)+SA+SO(IfK) 

W( I ) -Cl (K ) + DGK{ I ) +SYAw + CA*XZ(K) +SA+YZ IK ) 

220 G(I>KY+1/L)«G(I>KY-1>L )+<V(I)*( l.-H( I)+BZ(I )*HZ( I) )“U( I)*FXY(I )-W( 
1I)*HZ(I ) )/(FYY(I)*Bl(KY)} 

I-IO 

IF ( lO.NE. ITEKK) ) GO TO 230 
E»G(I3 jKY,L)-G(I0>KY#L) 

NO-NO+1 

E0(N0)»E0(N0)+P3*(E~E0(N0)) 

N»N0 

230 IF (I.LE.Il) GO TO 270 
I-I-l 
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E«0« 

IF ( IV( Ij.K) .NE.l) 60 TO 260 
ZZ-Z(K)-TYAW»(XC(K)+S1+A0( D+AOt I) ) 

240 IF (ZZ.GE.ZO(N-l) ) GO TO 250 
N»N-1 
GO TO 240 

250 RV»(ZZ-Z0(N-1) )/(Z0(N)-Z0(N-D) 
E»RV*E0(N)+(1 .~RV)*E0(N-1) 

260 M-NX+2-I 

G<I>KY + liL)»G(M#KY-ljU-E 
G(M,KY+1,L)-G(I>KY-1>L)+E 
GK2(M#KY)»GK1 (M>KY) 

GK1(M,KY)»G(M>KY#L) 

G(M>KY#U«G( I#KY»L) + E 
GO TO 230 
270 CONTINUE 

DO 280 I»2fNX 

280 G( I# J1~1>L ) -( l.-BETY/BSCALl)*G( I#J1>L) 
DO 290 J«lfMY 

6(I1-1> J»L )=( l.-BETX/6SCAL)*G(Il# J>L) 
290 6(I2+1#J>L)-(1.-BETX/BSCAL)*G(I2#J>L) 
RETURN 
END 
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